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00 ! Abstract 

Periodic normal forms for the codim 2 bifurcations of limit cycles up to a 3-dimensional center 
manifold in generic autonomous ODEs and computational formulas for their coefficients are derived. 
^0 ' The formulas are independent of the dimension of the phase space and involve solutions of certain 

, boundary- value problems on the interval [0,r], where T is the period of the critical cycle, as well as 

multilinear functions from the Taylor expansion of the right-hand sides near the cycle. The formulas 
allow us to distinguish between various bifurcation scenarios near codim 2 bifurcations. Our formulation 
, makes it possible to use robust numerical boundary-value algorithms based on orthogonal collocation, 

rather than shooting techniques, which greatly expands its applicability. The actual implementation is 
described in detail with numerical examples. 
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. 1 Introduction 

, Isolated periodic orbits (liniit cycles) of smooth differential equations 

(1) ii = fiu,p), u e M", p e R", 

play an important role in applications. In generic systems of the form ( 1 ) depending on one control parameter 
(i.e. with m = 1) a hyperbolic limit cycle exists for an open interval of parameter values p. At a boundary of 
■ such an interval, the limit cycle may become non- hyperbolic, so that either a cycle limit point (saddle-node) , 

, or a period- doubling (flip), or a torus [Neimark- Backer) bifurcation occurs. In two-parameter generic systems 

' (1) (i.e. with m = 2) these bifurcations happen at certain curves in the parameter plane. These curves of 

codim 1 bifurcations can meet tangentially or intersect transversally at some codim 2 points characterized 
by a double degeneracy of the limit cycle, which play the role of organizing centers for local dynamics, i.e. 
near the critical cycle and for nearby parameter values. In some cases, such codim 2 bifurcations imply the 
appearance of "chaotic motions" . 

The codim 2 bifurcations of limit cycles in generic systems (1) are well understood with the help of 
the corresponding Poincare maps and their normal forms (see for example [501 131 [Ml [IB])- However, 
applications of these results to the analysis of concrete systems (1) are exceptional, since they require 
accurate higher-order derivatives of the Poincare map which are hardly available numerically [HI [THl [HI [35] . 

We note that there exists software, e.g. CAPD [T], TIDES [H [3] that allows to compute up to any 
precision level the solution of an ODE using a Taylor series method in a variable stepsize - variable order 
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formulation. It can also compute, up to any order, the partial derivatives of the solution with respect to the 
initial conditions. When applied to compute a periodic orbit by a shooting method, this will also provide 
the derivatives of the Poincare map. 

Though this is a valuable approach in the case of a single periodic orbit, it is neither practical nor efficient 
in a continuation context since the periodic orbits are not obtained as the solution to a set of equations. Also, 
the shooting method does not have the high order convergence properties of the method of approximation by 
piecewise polynomials with collocation in the Gauss points that is routinely used in standard software such 
as AUTO [M], CONTENT [31], and MATCONT [1T1[T0]. Moreover, the number of derivatives of the Poincare 
map to be computed is 0{n'^) if derivatives up to order k are needed (in several cases fc = 5). Even for 
moderate values of n this involves a great deal of unnecessary work since in our situation the normal form 
itself is known in advance and we need only compute its coefficients. We will show that this can be done 
without computing the derivatives of the Poincare map. 

Indeed, recently an alternative numerical method to analyse codim 1 limit cycle bifurcations has been 
developed and implemented in [3D]. It is based on the periodic normalization proposed in [55] [531 HI] 
and completely avoids the numerical computation of Poincare maps and their derivatives. Instead, the 
computation of the normal form coeffcients is reduced to solving certain linear boundary value problems 
(BVP), where only the partial derivatives of the RHS of (1) are used [TTl dSl [2H] ■ In our implementation in 
MATCONT, we discretize these BVPs by orthogonal collocation with piecewise-polynomial functions. 

In the present paper, we apply the approach developed in [30] to codim 2 bifurcations of limit cycles. 
It should be noted that already in [^ normal forms for some codim 2 bifurcations of cycles in (1) were 
derived, while [23) contains the periodic normal forms for many codim 2 bifurcations of cycles, as well as 
a general normalization technique applicable at any codimension. However, in neither of these publications 
explicit formulas for the normal form coefficients were given. The derivation of such formulas is the primary 
contribution of this paper. 

The paper is organized as follows. In Section 2 we fix notation and list the periodic normal forms for 
codim 2 bifurcations of limit cycles. Then we derive explicit formulas to compute the critical normal form 
coefficients for these bifurcations, which we order by the dimension Uc of the cycle center manifold (i.e. the 
total number of critical multipliers with |^| = 1). We restrict to the cases ric = 2 and 3. The formulas are 
independent of the dimension of the phase space and involve solutions of certain BVPs on the interval [0, T], 
where T is the period of the critical cycle, as well as multilinear functions from the Taylor expansion of the 
right-hand sides of (1) near the cycle. A derivation of the critical periodic normal forms based on [23] is 
given in Appendix A, while their relationships with the normal forms of the Poincare maps are discussed in 
Appendix B 

2 Periodic normal forms on the center manifold 

Write ( 1 ) at the critical parameter values as 



and suppose that there is a limit cycle P corresponding to a periodic solution uq (t) = uo (i + T) , where T > 
is its (minimal) period. Develop F{uQ{t) + v) into the Taylor series 



(2) 



u = F{u), 



F{uo{t)+v) 



F{uo{t))+ 

A{t)v + -B{t; V, v) + -C{t; v, v, v) + 
^D{t;v,v,v,v) + ^E{t;v,v,v,v,v)+0{\\v\f), 



(3) 



where 



A{t)v = Fu{uo{t))v, B{t;Vi,V2) = Fuu{uoit))[vi,V2], C{t:Vi,V2,V3) = i^„„„(Mo(i))[«l, "2, ws], 
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etc. The multilinear forms A, B, C, D, and E are periodic in t with period T but this dependence will often 
not be indicated explicitly. 

Consider the initial- value problem for the fundamental matrix solution Y{t), namely, 



where /„ is the n x n identity matrix. The eigenvalues of the monodromy matrix M = Y{T) are called 
(Floquet) multipliers of the limit cycle. The multipliers with |^| = 1 are called critical. There is always a 
"trivial" critical multiplier fin — 1. We denote the total number of critical multipliers by ric and assume that 
the limit cycle is non- hyperbolic, i.e. ric > 1. In this case, there exists an invariant ric-dimensional critical 
center manifold W^'^(r) C R" near I^. 

It is well known [3l [29] , that in generic two-parameter systems ( 1 ) only eleven codim 2 bifurcations 
occur. To describe the normal forms of (2) on the critical center manifold W^'^(r) for these codim 2 cases, we 
parameterize VF^(r) near F by (r, ^), where r G [0, kT] for k G {1, 2, 3, 4} and ^ is a real or complex vector, 
depending on the bifurcation. It follows from [33] that it is possible to select the t- and ^-coordinates so 
that the restriction of (2) to the corresponding critical center manifold VF'^(r) with ric = 2 or ric = 3 will 
take one of the following periodic normal forms (for derivation see Appendix A ) . In Section 3 we will see 
that in two cases, namely the cusp of cycles bifurcation and the fold-flip bifurcation, a further simplification 
is possible in the normal form, and more specifically in the transformation of time. 

2.1 Bifurcations with a 2D center manifold 

We list here the critical periodic normal forms with Uc = 2 and briefly describe bifurcations in their generic 
unfoldings (see [21 [23 for more details). 

2.1.1 Cusp of cycles bifurcation 

The cycle has a cusp of cycles bifurcation (CPC) if the eigenvalue = = 1 of Y{T) corresponds to a two- 
dimensional Jordan block and the monodromy matrix has no other critical multipliers. The two-dimensional 
periodic normal form at the CPC bifurcation is 



where t e [0, T], ^ is a real coordinate on W'^{T) that is transverse to F, ai,a2,c £ M and the dots denote the 
0(^^)-terms, which are T-periodic in t. If c 7^ 0, the limit cycle F is triple. In generic two-parameter systems 
(1), three hyperbolic limit cycles exist in a cuspidal wedge approaching the codim 2 point and delimited by 
two bifurcations curves, where two cycles collide and disappear via the saddle-node bifurcation. 

2.1.2 Generalized period-doubling bifurcation 

The cycle has a generalized period- doubling bifurcation (GPD) if the trivial eigenvalue /i„ = 1 of the monodromy 
matrix Y{T) is simple and there is only another critical simple eigenvalue /ii = —1. The two-dimensional 
periodic normal form at the GPD bifurcation is 



(4) 



— =Ait)Y, Y{0)^I, 




(6) 




^This manifold should not be confused with the (ric — l)-dimensional center manifold of the corresponding Poincare map. 
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2.2 Bifurcations with a 3D center manifold 



where t S [0, 2T], ^ is a real coordinate on W^'^(r) that is transverse to F, ai, 02, e G M and the dots denote 
the 0(^^)-terms, which are 2T-periodic in r. li e ^ 0, at most two period doubled limit cycles can bifurcate 
from the critical limit cycle F. In generic two-parameter systems (1), the GPD-point in the period-doubling 
bifurcation curve separates its sub- and super-critical branches and is the origin of a unique saddle-node 
bifurcation curve, where two period doubled cycles collide and disappear. 



2.2 Bifurcations with a 3D center manifold 

We now list the critical periodic normal forms with — 3 and briefly describe bifurcations in their generic 
unfoldings (see [SllSn])- In all cases, "chaotic motions" are possible. 



2.2.1 Chenciner bifurcation 



The cycle has a Chenciner bifurcation (CH) if the trivial critical eigenvalue /x„ = 1 is simple and there are only 



two more critical simple multipliers /ii^2 = e with 9 7^ -j-, for j = 1,2,3,4,5,6. The three-dimensional 
periodic normal form at the CH bifurcation can be written as 



(7) 



(dr 
'dt 

dt 



:l + ai|eP + a2|ei' + 
ioj( + ic^\^\^+em' 



where r e [0,r], uj = 9/T, ^ is a complex coordinate on W^^(F) transverse to F, ai,a2,c G M, e G C and 
the dots denote the 0(||^^||)-terms, which are T-periodic in r. In generic two-parameter systems (1), at the 
CH-point the Neimark-S acker bifurcation changes its criticality (i.e. the bifurcating invariant torus changes 
its stability) . A complicated bifurcation set responsible for "collision" and destruction of two tori of opposite 
stability is rooted at this codim 2 point. 



2.2.2 Strong resonance 1:1 bifurcation 

The cycle has a strong resonance 1 : 1 bifurcation (Rl) if the trivial critical eigenvalue /i„ — 1 corresponds 
to a three-dimensional Jordan block. The three-dimensional periodic normal form at the Rl bifurcation is 



(8) 



^ = 6 + 66 + ... 

dt 



where r € [0,r], (616) a-re real coordinates on W'^{T) transverse to F, a,a,b e M and the dots denote the 
0(||^^||)-terms, which are T-periodic in r. In generic two-parameter systems (1), in the Rl-point is located 
on the saddle-node of cycles curve. At this point, a torus bifurcation curve is rooted together with global 
homoclinic bifurcation curves, along which the stable and the unstable invariant manifolds of a saddle cycle 
are tangent. The intersection of the invariant manifolds generates a Poincare homoclinic structure with the 
associated periodic and "chaotic motions" . 



2.2.3 Strong resonance 1:2 bifurcation 

The cycle has a strong resonance 1 : 2 bifurcation (R2) if the trivial critical eigenvalue /i„ = 1 is simple and 
there are only two more critical multipliers ^1 = /K2 = — 1 corresponding to a two-dimensional Jordan block. 
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The three-dimensional periodic normal form at the R2 bifurcation is 



(9) 



dr 
lit 

dt 



= l + a^l + ...., 

= 6 + aC?C2 + 



where r G [0,2T], (^1,^2) are real coordinates on W^T) transverse to F, a, a, 6 e M and the dots denote 
the 0(||^''j|)-tcrms, which are 2T-periodic in t. In generic two-parameter systems (1), the R2-point is the 
end-point of a torus bifurcation curve. The period-doubling bifurcation curve passes through this point, and 
(depending on the normal form coeffcients) a torus bifurcation curve of the period doubled limit cycle can 
originate there. As in the Rl-case, global bifurcation curves related to homoclinic tangencies can be present. 



2.2.4 Strong resonance 1:3 bifurcation 

The cycle has a strong resonance 1 : 3 bifurcation (R3) if the trivial critical eigenvalue = 1 is simple and 
there are only two more critical simple multipliers /ii,2 = e^'^. The three-dimensional periodic normal 
form at the R3 bifurcation can be written as 

^ = 1 + aiieP + + "3? + . • • , 

where r G [0, 3T], ^ is a complex coordinate on WIT) transverse to T, a\ G M, 02, as, 6, c e C with as = 0:2 
and the dots denote the 0(|]^'*||)-tcrms, which are 3T-periodic in t. In generic two-parameter systems (1), 
near the R3-point a homoclinic Poincare structure of the 3T-periodic limit cycle destroys the torus that is 
born at the Neimark-S acker bifurcation curve passing through this point. Homoclinic tangencies are rooted 
there. 



(10) 



2.2.5 Strong resonance 1:4 bifurcation 

The cycle has a strong resonance 1 : 4 bifurcation (R4) if the trivial critical eigenvalue /x„ = 1 is simple and 
there are only two more critical simple multipliers /ii,2 = e^'^ . The three-dimensional periodic normal form 
at the R4 bifurcation can be written as 

l + ai|e|2 + a2e^ + a3e^ + ..., 

c^\^\^+de + ..., 

where r e [0, 4T], ^ is a complex coordinate on W^T) transverse to F, ai € K, a2, as, c,d G C with as — 
and the dots denote the 0(||^^||)-tcrms, which are 4r-periodic in r. In generic two-parameter systems (1), 
at the R4-point there can be eight different situations, depending upon the values taken by the parameter c 
and d. In the simplest case a homoclinic structure associated to a 4T-periodic cycle destroys an invariant 
torus that is born at the Neimark-S acker bifurcation curve that passes trough this point. 

2.2.6 Fold-flip bifurcation 

The cycle has a /o/d flip bifurcation (LPPD) if the trivial critical eigenvalue /ii = /i„ = 1 is double non 
semi-simple and there is only one more critical multiplier /i2 = —1. The three-dimensional periodic normal 



(11) 
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form at the LPPD bifurcation is 



' dT 



1 - Cl + "20^1 + "02^2 + "30^1 + ai2ClC2 + ■ • ■ . 



(12) 



< 



da 



a20Cl + "02^2 + ^30^1 + 012^1^2 + ■ ■ • 



dt 
. dt 



where t G [0,2T], (6j6) are real coordinates on M^'^(r) transverse to F, all the coefficients are real and 
the dots denote the 0(||^''||)-terms, which are 2T-periodic in t. In generic two-parameter systems (1), the 
period-doubling and saddle-node of cycles bifurcation curves are tangent at the LPPD-point, where (depending 
on the normal form coefficients) a Neimark-Sacker bifurcation curve of the 2r-periodic cycle can be rooted. 
Global bifurcations of heteroclinic structures and invariant tori are also possible. 

3 Computation of critical coefficients 

In view of the above, we can assume that a parametrization of the center manifold ly^(r) has been selected 
so that the restriction of (2) to this manifold has one of the normal forms given in Section 2, We then 
apply the so-called homological equation approach [6]: the Taylor expansions of T-, 2T-, 3T- or 4T-periodic 
unknown functions involved in these parametrizations can be found by solving appropriate BVPs on [0, T] 
so that (2) restricted to ly^(r) has the corresponding periodic normal form. The coefficients of the normal 
forms arise from the solvability conditions for the BVPs as integrals of scalar products over [0,T], involving 
nonlinear terms of (2) near the periodic solution uq, as well as the critical eigenfunctions and higher order 
expansion terms of the center manifold. The Taylor expansion coefficient functions are usually unique up to 
the addition of a multiple of a known eigenfunction. This has to be fixed by adding an integral condition. 
Among other things this leads to the fact that normal form coefficients are not unique but implications 
for the underlying dynamical systems are independent of this. We also remark that the solvability of all 
the equations up to the maximal order of the normal form has to be checked. Finally, we note that the 
coefficients related to the transformation of time will only be computed when needed in the computation of 
the critical coefficients in the normal form for the state variables. 

3.1 Cusp of cycles bifurcation 

The two-dimensional critical center manifold Ty^(r) at the CPC bifurcation can be parametrized locally by 
(t, as 



where uq and all hj are functions of r, with hj{T) = hj{0), for j — 2, 3, while the generalized eigenfunction 



The function v exists due to Lemma 2 of |23| . Let ip* be a nontrivial solution of the adjoint eigenvalue 
problem 




V is given by 



(15) 



V - A{t)v - F{uo) = 0, 
v{T)-v{0) = 0, 

{v,Fiuo))dT = 0. 




(16) 




0, TG [o,r] 

0, 
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and the generalized adjoint eigenfunction v* a solution of 
(17) 

which is now defined up to the addition of a multiple of (p*. Note that the first equation of (15 ) implies 



i!* + A^{t)v* + ip* = 0, TG [o,r], 

v*{T)-v*{0) = 0, 



(18) r{<fi*,F{uo))dT = 

Jo 



'0 

for ip* satisfying (16). Indeed, 

rT i-T 



^ {v*,Fiuo))dT^ {^*^i^±-A{T)^v)dT = ~ {(^^+A^{r)^^*,v)dT = Q 
due to (16). 

Moreover, due to spectral assumptions at the CPC-point, we can also assume 

(19) / {^%v)dT = l. 

Jo 

Notice that this assumption gives us another normalization for free, since because of (15) and (17) we have 

(20) 1^ {v*, F{uo))dT = £ {v*, (^^ - A{T)j v)dT 



= - j\(^^+A^(T)^v\v)dT = j\<p\v)dT=l, 



i.e. wc have normalized the eigenfunction of the adjoint problem with the generalized one of the original 
problem and the generalized eigenfunction of the adjoint problem with the eigenfunction of the original 
problem. So ip* is the unique solution of the BVP 

r ^* + A^{t)^* = 0, re [o,r], 

(21) <^ ^*{T)-ip*{Q) = 0, 

i J^{<p*,v)dT-l = 0. 

Now, wc still nccid an integral condition for the adjoint generalized eigenfunction v* . In all cases, for 
the computation of an adjoint generalized eigenfunction we will require the inproduct with an original 
eigenfunction to be zero. Here, the inproduct with v is appropriate. Therefore we obtain 

r V* +A'^{t)v* +(p* = 0, re [0,T], 

(22) I v*{T)-v*{0) = 0, 

I Io{v*,v)dT = 0. 

Now, we substitute (13) into (2), using (3), (5) and (14), as well as 

du du d^ du dr 
lit ~ d^di ^ drlt' 

This gives 

iio+^iv - uq) + S,^ (ctiuo - V + 7:^2) + [c(2Uo + ctiv - i/i2 + ^/i3 + cv) + 0{£,^) 



= F{uo) + ^A{t)v + iA{T)h2 + B{t; v, v)) + If (A(t)/i3 + 3B(t; h2,v) + C(r; v, v, v)) + 0{f), 
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3.1 Cusp of cycles bifurcation 



where dots denote the derivatives with respect to r. 
Collecting the ^°-terms we get the identity 



iio = F{uo), 

since uq is the periodic solution of (2). 

The -terms provide another identity, namely 

i) — iiQ — A{t)v, 

due to (15). 

From collecting the ^^-terms we obtain an equation for /12 

(23) /i2 - A{T)h2 = B{t; v, v) +2v- 2aiUo. 

Now, we project the left-hand side of this equation on the adjoint null-eigenfunction and obtain 

^{^*, (^^ - A(t)) h^) dT^-j\iy^^+ A^{t)^ dr = 

due to (16). We can use this result to impose the so-called Fredholm solvability condition, i.e. also project 
the right-hand side of (23) on Lp* which then has to be equal to 

{if* , B(r; u, v)+2v~ 2aiUa) dr = {tp* , B{t; v, v) + 2A{t)v) dr ^ 0. 

Jo 

Notice that this condition is actually trivially satisfied, due to the fact that we are at a cusp of cycles point, 
so that the second order normal form coefficient 

1 

- {ip*,B{T;v,v) + 2A{t)v) dr 
2 Jo 

(see |30j) vanishes. Hence equation (23) is solvable, independent of the value of ai. So, for any value of ai 
we get an equation for h2 to be solved in the space of vector-functions on [0,T] satisfying /i2(T) = /i2(0). 
Notice that if /12 satisfies (23), /12 + £F{uo) also satisfies (23), due to the fact that F{uo) is in the kernel 
of the operator ^ — A{t) and to the linearity of this operator. Now, the orthogonality condition with v* 
determines the value of e such that we can define /12 as the unique solution of 

r h2- A{T)h2- B{t;v,v) -2Av-2F{uo)+2aiF{uo) = 0, Te[0,T], 

(24) I /i2(r)-/i2(0) = 0, 

[ J^{v*,h2)dT = 0. 

Collecting the ^'^-terms we finally obtain an equation in /13 which allows us to determine the normal form 
coefficient c of (5) 

/13 — A(t)/i3 = —6a2UQ — Qaiv + 3/i2 — 6cu + •iB{T; /12, v) + C{t; v, v, v). 

As before, the null-eigenfunction of the adjoint operator — ^ — ^'^i''') is V*- Thus, the Fredholm solvability 
condition implies that 

T 

(ip* , —602^0 — Gctiv + 3/i2 — 6cv + 3B{t; h2,v) + C{t; v, v, v)) dr = 0. 



Using (15), (19) and (18), we get the expression 
6 Jo 



(25) c^l I {ip*,-6aiA{T)v + 3A(t)/i2 + 3B(r; u, v) + 6A{t)v + 3B{t; h2,v) + C{t; v, v, v)) dr, 
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where v and (p* are defined by (15) and (21), while /i2 satisfies (24). 

In what follows we will prove that the choice of ai does not influence the value of the critical normal 
form coefficient c. Two solutions /12 differ by h'^^ — h^2^ ~ —2(0;^^ — a^^^)v, from which 



-(2)-c(i) = J I -6(af ^ -a^'VM^'-6A(T)(af^ -ai'V-6(af^ -a^'^)B(T;?;,?;)) dr 
Jo 



T 



(af^-aj^^) / {if\-2A{T)v- B{t-v,v)) dT 







T 



= {af^-at^)l {i{>\-2A{T)v + 2A{T)v) dT 
= 0, 







since b vanishes. So, for simplicity we take ai = which further simplifies normal form (5). 

Therefore, the critical coefficient c in the periodic normal form (5) has been computed. The bifurcation 
is nondegenerate if c 7^ 0. 

3.2 Generalized period-doubling bifurcation 

The two-dimensional critical center manifold Vl^'^(r) at the GPD bifurcation can be parametrized locally by 
(t, as 

(26) u^u^ + ^v^ h{t,0. r e [o,2r], e e 

where the function H satisfies H{2T,^) ~ H{0,S^). It has the Taylor expansion 

(27) H{t, = \h2e + \h:ie + Y^h^e + ^^^^5^' + o{e), 
with hj{2T) = ft.j(0), with 

v-A{t)v = 0, T e [0,T], 

(28) { v{T) + viO) = 0, 



and 



J^{v,v)dT-l = 0, 



v{t + T) = ~v{t) for T £ [0, T]. 



The function v exists due to Lemma 5 of [23j . 

The functions hi, i = 1 . . .5, can be found by solving appropriate BVPs, assuming that (2) restricted to 
W^(r) has the periodic normal form (6). From (26) and (27) it follows that hi{T + T) = hi{T) for i even 
and hi{T + T) = —hi^r) for i odd, for r G [0, T]. Indeed, since we are at a generalized period-doubling point 
w(t,0 = w(r + T,-e), so 

i i 

and thus 

from which the stated follows. This makes it possible to restrict our considerations to the interval [0, T] 
instead of [0,2T]. 

The coefficients ai, 02 and e arise from the solvability conditions for the BVPs as integrals of scalar 
products over the interval [0,T]. Specifically, these scalar products involve among other things the terms 
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3.2 Generalized period-doubling bifurcation 



up to the fifth order of (1) near the periodic solution uq, the eigenfunction v, the adjoint eigenfunction (p* 
satisfying 

r r+A^{T)<f* = 0, re[0,T], 

(29) I ^*{T)-ip*{0) = 0, 

[ J^{ip*,F{uo)}dT-l = 0, 

and a similar adjoint eigenfunction v* satisfying 

r v*+A^{t)v* = 0, re [0,T], 

(30) I v*{T) + v*{0) = 0, 

i {v*,v)dT-l = 0. 

To derive the normal form coefficient, we proceed as in Section 3.1j namely, we substitute (26) into (2) 
and use (3), as well as (6) and (27). 

Collecting the ^"^-terms in the resulting equation gives us the identity 

uo = F{uo), 

since uq is the T-periodic solution of (2). 
The -terms provide another identity 

V = A{t)v, 

due to (28) and (3.2). 

By collecting the ^^-terms, we obtain the equation for /12, 

(31) h2 - A(r)/i2 = -B(t; v, v) - 2aiuo, 

to be solved in the space of functions satisfying h2{T) = /i2(0). In this space, the differential operator 
^ — A{t) is singular with null- function Uq. Thus, the following Predholm solvability condition is involved 



1 



which leads to the expression 



T 

{tp* , B[t;v,v) — 2aiUo) dr = 0, 



T 



(32) "'"2 70 

where V and cp* arc defined by (28) and (29), respectively. 

With ai defined in this way, let /12 be a solution of (31) in the space of functions satisfying /i2(0) = h^iT). 
Notice also here that if /12 is a solution of (31), then also /i2 + £iF{uo) satisfies (31), since F{uq) is in the 
kernel of the operator ^ — A{t). In order to obtain a unique solution (without projection on the null 
eigenspace) we impose the following orthogonality condition which determines the value of ei 



/ (^*,/i2) rfr = 0. 
Jo 



Thus /i2 is the unique solution of the BVP 

r h2- A{T)h2- B{T;v,v) + 2aiF{uQ) = 0, rG[0,r], 

(33) <^ /i2(T) - /i2(0) = 0, 
i !^{^*M)dT = 0. 

Collecting the ^^-terms, we get the equation for /13, 

(34) hi - A(r)/i3 = C(r; v, v, v) + 3S(r; v, /12) - Qa^v, 
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to be solved in the space of functions satisfying h-^iT) — —h-^iQ). In this space the differential operator 
^ — A{t) has a one-dimensional nullspace, spanned by w, and (34) is solvable only if the RHS of this 
equation lies in the reachable space of that operator. Using (28), we can rewrite the right-hand side as 

C(t; w, w, v) + SB{t; v, /12) — 6aiA{T)v. 

But the Fredholm solvability condition 

(35) / {v*, C{t; V, V, v) + 3B{t; V, h2) - 6ai A{t)v) dT = 



is trivially satisfied due to the fact that we are in a generalized period-doubling point and so the cubic 
coefficient of the normal form 

1 '•^ 



c=- / (v* ,C(t;v,v,v) + 3B(T;v,h2) — 6aiA{T)v) dr 
6 Jo 

(for it's definition see [30]) vanishes. Since the RHS of (34) is in the range space of the operator ^ — A{t), 
we can solve the equation in order to find /13 as the unique solution of the EVP 

( hs- A{T)h3~C{T;v,v,v) -3B{T;v,h2) + 6aiA{T)v ^ 0, t€[0,T], 
(36) <^ /i3(r) + /i3(0) = 0, 

[ !^{v*M)dT = 0. 

By collecting the ^''-terms, we get the equation for /14, 
/14 - A(t)/i4 = £)(r; u, u, u, v) + 6C(t; w, w, /12) -|- 3B(t; /12, '^2) 

+ 4B{t] V, /13) - 12q!i/i2 - 24a2ito, 

to be solved in the space of functions satisfying h4{T) = ^,4(0). Formulation of the Fredholm solvability 
condition 

T 

{iP*,D{t;v,v,v,v) +6C{T;v,v,h2) -(- 3B(r; ftg, ^2) +4B(t;u, /13) - 12ai/i2 - 24a2^to) dr ^ 



gives us an equation for the parameter a2 



1 '•^ 



Q!2 = — / {ip*,D{T; v,v,v,v) + 6C(r; v, v, ^2) + 3B{t; /i2, /12) + 4i3(T; v, ^13) - 12ai/i2) rfr 
24 Jo 

which can be simplified considering (31 ) into 
1 '-^ 



012 = — / {ip*,D{t;v,v,v,v) +6C{T;v,v,h2) + 3B{T;h2,h2)+ 
^4 Jo 

4B{T;v,h3)-12ai{A{T)h2+B{T;v,v))) dr + aj, 

where ai is given by (32), and /12, /13, w and (yS* are the solutions of the BVPs (33), (36), (28) and (29), 
respectively. 

Using this value of a2 we can find ft.4 by 



(37) 



/14 — A{T)h4 — D{t; V, V, V, v) — 6C(r; v, v, ft.2) — 3i?(T; /12, /12) — 

4S(t; w, /i3) + l2ai{A{T)h2 + B{t; v, v) - 2qiF(mo)) + 24a2F(uo) = 0, r G [0, T], 

/l4(T) - /l4(0) = 0, 

/o^(¥.*,MrfT = 0. 
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Finally, by collecting the ^^-terms, we get the equation for h^, 

/15 — A(t)/i5 — E{t; V, V, V, V, v) + 10D{t; v, v, v, /12) + 15C(t; v, /12, ^2) 

+ 10C(t; V, V, hs) + 10B{t; /i2, hs) + 5B{t; v, /14) - 12002" - 20ai/i3 - 120eu, 

which has to be solved in the space of functions satisfying h^{T) — —h^{0). Since the operator ^ — ^(''") 
has a one-dimensional null-space, we can write 

i-T 

{v* ,E{t; V, V, V, V, v) + 10D{t; v, v, v, /12) + 15C(t; v, /12, ^12) + 10C(t; v, v, /13) -|- 



10B{t; /i2, /la) + 5B(r; v, hi) - 120a2i' - 20ai/i3 - 120eu) dr = 0, 

which makes it possible to compute the parameter e of the normal form (6). Using the normalization of 
(30), (36) and (35) gives 

1 

(38) e=— / {v*,E{t;v,v,v,v,v) + 10D{T;v,v,v,h2) + l5C{T;v,h2,h2) + 

10C(t; V, V, hs) + 10B{t; /12, /13) -I- 5B{t; v, /14) - 120a2^(T)w - 2QaiA{T)h^) dr. 
A check that this quantity doesn't vanish, guarantees us that the codim-2 bifurcation is non-degenerate. 

3.3 Chenciner bifurcation 

The three-dimensional critical center manifold W^'^(r) at the CH bifurcation can be parametrized locally by 
(t, C) as 

(39) u = u^ + iv + lv{T)+H{T,^,0, Te[0,T], ^ec, 

where the real function H satisfies H{T,£^,^) = i/(0,^,^), and has the Taylor expansion 

+ Ihsof + lh2iei+ ^/»i2^f + IhosC" 

+ ^hMr)e + lh,,{T)e^+ h22ir)e^ + ^/ii3(r)ef + ^/i04(r)f 

+ Y^^soM^' + ^^4i(r)e^e + ^h^2ir)e^' + ^/^23(r)C^f -I- ^h^^rM'' 

(40) + ^ho,{r)^+Oim, 

with hij{T) = hij(0) and hij = hji so that ha is real, while v and its conjugate v are defined as 

( v{t) - A{t)v + iuj V = 0, re[0,r], 

(41) } v{T)~viO) ^ 0, 

[ J^{v,v)dT-l = 0. 

These functions exist due to Lemma 2 of p3]. 

If we assume that (2) restricted to Ty^(r) has the periodic normal form (7), as in the previous cases, we 
can find the functions hij{T) by solving appropriate BVPs. 

First we introduce the two needed adjoint eigenf unctions. The first one, namely ip*, satisfies (29), and 
the second one, namely v* , satisfies 

i* {t) + {t)v* + iuj V* = 0, re[0,T], 

(42) { v*iT)^v*iO) = 0, 

{v*,v)dT-l = 0. 
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As usual, we substitute (39) into (2), use (3), (7), and (40), as well as the homological equation 

du du du d^ du dr 

and collect the corresponding terms in order to find the needed coefScients of (7). 
The ^-independent and the linear terms give rise to the usual identities 

Uq = F{uq), V — A{t)v + iijjv = 0, V — A(t)v — itov = 0. 

Collecting the coefficients of the or ^^-terms leads to the equation 

h20 - A(r)/i20 + 2iujh2o = B{t; v, v) 

or its complex-conjugate. This equation has a unique solution /120 satisfying /i2o(T) = h2o{0), since due to 
the spectral assumptions e^**^^ is not a multiplier of the critical cycle. Thus, /120 can be found by solving 

(.^^ i h2o-A{T)h2o + 2ituh2o-B{T:v,v) = 0, Te[0,T], 

^ ' 1 /i20(T) - /120(0) = 0. 

By collecting the ^^-terms we obtain an equation for hn, namely 

/ill - A{T)hn = B{t; v, v) - aiuo, 

to be solved in the space of the functions satisfying hn{T) = /in (0) . In this space the operator ^ — A{t) has 
a range space with codimcnsion one. As before, the null-eigenfun(;tion of the adjoint operator — ^ — A'^{t) 
is (f*, given by (29), and thus because of the Fredholm solvability condition, we can easily obtain the needed 
value for ai 

(44) ai= [ {ip*,B{T;v,v))dT. 

Jo 

With ai defined in this way, let hn be the unique solution of the BVP 

r hn- A{T)hii- B{T;v,v) + aiUo = 0, re[0,T], 

(45) I hn{T)-hn{0) = 0, 

I j^{^*,hn)dT = 0. 

The coefficient in front of the third order terms in (7) is purly imaginary since the first Lyapunov 
coefficient vanishes at a Chenciner point. We are now ready to compute this coefficient. In fact, if we collect 
the ^^f- terms we obtain 

h2\ — A{T)h2i + iojh2\ = C{t; v, v, v) + 2B{t; v, hn) + B{t; v, /120) - 2icv - 2aiv, 

to be solved in the space of functions satisfying h2i{T) = ft.2i(0). In this space the operator ^ — A{t) + iui 
is singular, since e*"^ is a multiplier of the critical cycle. So we can impose the usual Fredholm solvability 
condition, taking (42) into account 

(46) [ {v*,C{T;v,v,v) + 2B{T;v,hn)+B{T;v,h2o)-2icv-2aiv)dT = 0. 
Jo 

From this equation we can find the value of the coefficient c of the normal form (7) 

(47) •==--/ {v\C{T;v,v,v) + 2B{T;v,hn)+B{T;v,h2o)-2axA{T)v)dT + axuj 

2 Jo 
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and, with c defined in this way, we can find ft,2i as the unique solution of the EVP 



(48) 



/i2i - ^(r)/i2i + icj/i2i - C{t: v, v, v) - 2B{t; v, hn) 

-B{T;v,h2o) + '2icv + 2ai{A{T)v-iuiv) = 0, re[0,T], 

/i2i(r)-/i2i(o) = 0, 

Jo{v*,h2i)dT = 0. 



Collecting the ^^-terms gives us an equation for /130 

/130 - A{T)h3o + Siujhso = C'{t; v, v, v) + 3B{t; v, /120), 

which has a unique solution ^,30 satisfying h3o{T) = /i3o(0), since e"^*"-^ is not a multiplier of the critical 
cycle by the spectral assumptions. Thus, h^o is the unique solution of the EVP 



(49) 



iiso- A{T)h3o + 3iwh3o-C{T;v,v,v) -3B{T;v,h2o) = 0, re[0,T], 

h3o{T) - h3o{0) = 0. 



By collecting the ^^^-terms we obtain an equation for /131 

A31 - A(t)/i3i + 2iw/i3i = D{t;v,v,v,v) + 3C{T;v,v,hn) + 3C{t;v,v, h2o) + 3B(t; /in, /120) 

+3B(t; V, /121) + B{t; v, /130) - 6ic/i2o - 3ai/i20 

which has a unique solution ^,31 satisfying /i3i(T) = /i3i(0), since e^*"^ is not a multiplier of the critical 
cycle by the spectral assumptions. Thus, /131 is the unique solution of the EVP 



(50) 



A31 - A(r)/i3i + 2iLoh3i - D{t: v, v, v, v) - 3C(r; v, v, hn) 
-3C(t; u, V, /i2o) - 3B(r; /in, /120) - 3B(r; v, /121) 
-B{t; V, hso) + 6ic/i2o + 3ai (A(t)/i2o - 2iu!h2o + B{t; v, v)) 

/l3l(T)-/l3l(0) 



0, re [0,T], 
0. 



Taking into account the |^|^-terms gives an equation for /122 



/122 - A(r)/i22 = D{t;v,v,v,v) + C{T;v,v,ho2) 

+4C(r; t;, /in) + C(t; u, /120) + 25(r; /in, /in) + 2B{t; v, /112) 
+B(t; /102, /120) + 2i?(r; w, /121) - Aaihn - 402^0, 

to be solved in the space of functions satisfying h22{T) = /i22(0). In this space the operator ^ — A{t) has 
a range space with codimension one which is orthogonal to ip*. So one Predholm solvability condition is 
involved, namely 



Jo 



{iP*,D{t;v,v,v,v) + C{T;v,v,ho2) + 4C{T;v,v,hii) + C{T;v,v,h2o) + 2B{t; /in, /in) 



+ 2B{t; V, /112) + B{t; /102, /120) + 2B{t; v, /121) - 4aihn - 4a2Uo)dT = 0, 
which allows us to compute the value of the coefficient a2 of our normal form 

1 

(51) a2 = - {(fi*,D{T;v,v,v,v) + C{T;v,v,ho2)+4:C{T;v,v,hn) 
4 Jo 

+ C(t; V, V, /i2o) + 2B{t; /in, /in) + 2B{t; v, /112) + B{t; /102, /120) 

+ 2B{t; V, /i2i) - Aai{A{T)hn + B{t; v, v)))dT + af . 
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Using this value for a2 we can find /122 as the unique solution of the BVP 



(52) 



/i22 - A(t)/i22 - D{t; V, V, V, v) - C{t; v, v, /102) 
-AC{t;v,v, /ill) - C{T;v,v,h2o) - '2'B{t; /in, /in) - 2B(r; /112) 

-B{t; ho2, /120) - 2-B(t; w, /121) 
+4q;i(A(t)/iii + B{t; v, v) - aiF{uo)) + 4a2F{uo) 

h2l{T) - /l22(0) 



0, TG [0,T], 
0, 



j^{ip*,h22)dT = 0. 



Finally, by collecting the ^^^^-terms we get an equation for /132 

h32 — A{T)h32 + iw/i32 = E{t; v, v, v, v, v) + D{t-, v, v, v, /102) + 6D(t; v, v, v, hn) 
+3D{t; V, V, V, /120) + 6C(t; v, hn.hn) + 3C(t; v, v, hi2) + 3C(t; w, /102, /120) 
+6C(t;'u,/iii,/i2o) + 6C(t; -y, -y, /121) + C(t; i), w, /i3o) + 3B(t; /112, /120) + 6-B(t; /in, /121) 

+35(r; w, /122) + S(t; /102, /130) + 25(r; w, /131) 
— 12et; — 6zc/i2i — 12a2V — 60:1/121 

that, since the operator is singular, allows us, using the first of (41) as well as the first and the last of (48) 
and (46), to compute the critical coefficient e of (7) imposing the Predholm solvability condition, obtaining 

1 /"^ 

E{t; V, V, V, V, v) + D{t; v, v, v, /102) + 6£>(r; v, v, v, hn) 



+3D{T]v,v,v,h2o) + 6C{T;v,hii,hii) + 3C(r; /112) 
+3C(t;w,/io2,/i2o) +QC{t;v, /iii,/i2o) + 6C{T;v,v,h2i) + C{T;v,v,h3o) 
+3B{t; hi2, h2o) + 6B{t; hnMi) + 3B(t; v, /122) + S(r; /102, /130) 
+2B{t- V, hsi) - 12a2A{T)v - 6ai(A(T)/i2i + 2B(r; v, /in) + C(r; v, v, v) 

(53) +B{t; V, /120) — 2aiAv))dT + iuja2 + icai — cx{iuj. 

As stated in Appendix B, we define the second Lyapunov coefficient as 

£2(0) =3ff(e). 

If this coefficient does not vanish, no more degeneracies happen at this codim 2 point. 

Since we have to check all equations up to the fifth order, we still have to look at the ^^-terms, the 
5^-terms and the ^^^-tcrms, which give respectively 



/140 — A(t)/140 + 4iw/l40 

/150 - A(r)/i5o + 5iu;/i5o = 



= D{t] V, V, V, v) + 6C(r; v, v, /120) + 3B{t; /120, /120) 
B{T;v,h3o), 

E{t; V, V, V, V, v) + 10D(t; v, v, v, /120) + 10C(t; v, v, /130) 
+15C(r; V, /120, /120) + 10-B(r; /120, /130) + 5-B(t; v, /140) 



and 



/141 — A(t)/i4i + 3iu)h4i = E{t; v, v, v, v, v) + 6D{t; v, v, v, /120) + 4Z)(r; v, v, v, hn) 

+4C(t; V, V, /130) + 8C(r; v, /120, hn) + 6C(t; w, w, /121) 
+C(t; w, /120, /120) + GB{t; /120, /121) + B{t; v, /140) 
+4i3(r; /131) + 4B(r; /130, /in) - 3qi/i3o - 12ic/i3o. 
No solvability conditions have to be satisfied. 

Since wc; arc; in a complex eigenvalues case, v is determined up to a factor 7, for which 7^7 = 1. Then 
V*, /120, /121, /i30i hsi are replaced by jv*, 7^/120, 7/121, 7'^/i30, 7^/i3i respectively, but ai,a2,c and e are not 
affected by this factor. 
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3.4 Strong resonance 1:1 bifurcation 

The three-dimensional critical center manifold Ty^(r) at the Rl bifurcation can be parameterized locally by 
(r, ^) as 

(54) u^uo+iiVi+^2V2+HiT,0, r e [0,T], e= (6,6) eK', 
where H satisfies H{T,^) = H{0,^) and has the Taylor expansion 

(55) h{t, 6 = h^ofi + /»ii66 + lho2ii + OiUf), 

where the functions /i20j ^ii and /102 are T-periodic in r, where vi and V2 are the generalized eigenfunctions 
associated with the trivial multiplier and defined as the unique solutions of the BVP 

r ~ A{t)vi - F{uo) = 0, re [0,T], 

(56) <^ vi{T)~vi{0) = 0, 

and 

( V2-A{t)v2+Vi = 0, Te[0,T], 

(57) <^ W2(T)-W2(0) = 0, 

1 {v2,Fiu„))dr = 0, 

respectively. The functions vi and V2 exist and are different due to Lemma 2 of [23] . Following our approach 
to find the value of the normal form constants, we define if* as a solution of the adjoint eigenfunction problem 
(16), vl as a solution of the adjoint generalized eigenfunction problem (17) and V2 as a solution of 

" v2*iT) + A^{t)v*2 +vt = 0, T e [o,r], 

v*2{T) ~ vm = 0. 

First, notice that the Fredholm solvability condition gives us immediately the following scalar products 

(58) / {ip*,Fiuo))dT= [ {ip*,vi)dT= [ {Fiuo),vl)dT^O. 
Jo Jo Jo 

Due to the spectral assumptions at the Rl point we are free to assume that 

(59) / {ip*,V2)dT^l. 

Jo 

Appending this condition to the eigenproblem, we can find the eigenfunctions ip* as the unique solutions of 
the BVP 

r r +AT(r)^* = 0, re [0,T], 

(60) <^ ^*(T)-^*(0) = 0, 

[ /„^(^*,t;2)dr-l = 0. 

As already mentioned in the cusp of cycles case, we will choose adjoint generalized eigenfunctions orthog- 
onal to an original eigenfunction. Therefore, vl and V2 are obtained as the solution of 

vi* + {t)vI - = 0, re[0,T], 

(61) { vl{T)-vm = 0, 

/o^(«i:'^2)dr = 0, 
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and 

r «2*(t) + A^{t)v*2 +vt = 0, r G [0,T], 
(62) I v*^{T)-vm = 0, 

I lo {V2,V2)dT = 0, 

respectively. Notice that, as in the cusp of cycles case, we have normalized in (59) the adjoint eigenfunction 
with the last generalized eigenfunction, which gives us in addition 



/ {vl,vi)dT= {v2,F{uo))dT = 1. 
Jo Jo 



As usual, to derive the value of the normal form coefficients wc substitute (54) into (2), wc use (3) as 
well as (8) and (55) and get different equalities for every degree of ^. Remark that in fact the solvability of 
all the equations up to the maximal order of the normal form has to be checked. We will pay extra attention 
to it in this section. 

By collecting the ^°-terms we get the identity 

Uo = F{uo). 

The linear terms provide two other identities, namely 

vi — A{t)vi — F{uo) = 0, V2- A.V2 +vi =0, 

cf. (56) and (57). 

By collecting the ^f-terms we find an equation for /120 

(63) /120 - A{T)h2o = -2q!Uo + 2*1 + B{t;vi,vi) - 2av2, 

to be solved in the space of periodic fimctions on [0,T]. In this space, the differential operator ^ — A(t) is 
singular with a range orthogonal to (p* . Thus a Fredholm solvability condition is involved, namely 



L 



T 

(if*, —2auo + 21)1 + B{t; vi,vi) — 2av2)dT = 0. 



'0 

The equations (58), (59), and (56) let us obtain the following value for a 

(64) «=J/ {v*,'^A{T)vi+B{T;vi,vi))dT. 

^ Jo 

Notice that in the RHS of (63) we have no freedom which could change the value of the coefficient a. This 
confirms the theoretically proved fact that the ^^-term of normal form (8) is resonant. Notice moreover that 

parameter a is undetermined, which gives us two degrees of freedom for ft,20- In fact, if ft,20 is a solution of 
(63), then also /120 = ^20 + £20-^(^0) + ^io'^i is a solution, due to the fact that F{uo) spans the nuUspace of 
the operator ^ — A{t) and that we can tune a as desirable: 

•^^20 .7 dh2o , II ( dvi \ dh2o ^, , 11 ■ 



(65) - A{T)h2o = -^- ^(^)/i2o + £20 - A{T)v^j = ^ - ^(^)'^2o + 4o«o. 
By collecting the ^1^2-terms we find an equation for h\\ 

(66) /ill - A{T)hii = B{t-, vi,V2) + 1)2 - h2o - bv2 - wi, 

to be solved in the space of T-periodic functions. As in the previous case, a solvability condition is involved 



I 



T 

{iP*,B{t;vi,V2) +V2- h2o - bv2 - vi)dT = 0. 
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Equation (59) as well as (57) and (58) let us rewrite this condition as 

b= [ {^*,B{T;vi,V2)+A{T)v2)dT- [ {ip* ,h2o)dT. 
Jo Jo 

Using (61), (63), (58) and (61) we can rewrite the second term of the right-hand side as 

j\(J^+A^{T)^ vl,h20)dT = j\vl, (^-1-+A{t)^ h20)dT 

T pT 

{vl,-2auo + 2vi + B{t;vi,vi) - 2av2)dT = - / {vl,2Avi + B{T;vi,vi)}dT 



Jo 

obtaining an equation for b which involves only the original and adjoint eigenfunctions 

(67) b= [ {^* ,B{t;vi,V2) + A{T)v2)dT + [ {vl,2Avi + B{T;vuVi))dT. 

Jo Jo 

Notice that the freedom that we have on /120 can not be used to change the value to coefficient b (and so 

the of the normal form (8) is resonant). Indeed, ^20 is defined up to a multiple of F(uo) and 

vi, but both vectors are orthogonal to if*, see the first two orthogonality conditions in (58). However the 
presence of /?2o in the RHS gives us three degrees of freedom for hn- In fact, if hn is a solution of (66), also 
hn = hn + e{iF{uQ) — £20^1 + ^20^2 is a solution, since 



dhi 



dt 



= - A{T)hii - eioF{uo) - £20^1 • 

Collecting the Q-terms gives us the following equation for ho2 

ho2 - A{T)ho2 = B{t, V2,V2) - 2hn, 

to be solved in the space of T-periodic functions. This equation should be solvable, so the RHS should lay 
in the reachable space of the operator ^ — A{t): 



L 



T 

{ip*,B{T,V2,V2)-2hn)dT = 0. 



10 

This condition can be satisfied by correctly tuning hn. In fact, £30 is not yet determined, so hu can have a 
projection on V2- Due to(59) V2 does not lay in the reachable space of the ^ — A{t) operator, and therefore 
we can impose that 



j {ip*,hn)dT =^ j {ip*,B{T,V2,V2))d7 



This last solvability condition determines £20 uniquely, and since £30 determines the value of a, see (63) 
and (65), also a is now uniquely determined. So the center manifold expansion (8) has now become unique. 
Note that in fact the value of a is not needed since, as shown in Appendix B , it does not affect the bifurcation 
scenario. Remark also that in order to compute the necessary coefficients a and b by equations (64) and (67), 
the second order expansion of the center manifold is not needed. Indeed, we have rewritten the formulas of 
the normal form coefficients in terms of the original and adjoint eigenfunctions. /120 or hn are not needed, 
therefore we don't write down the EVPs for their unique solutions. 
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3.5 Strong resonance 1:2 bifurcation 

The three-dimensional critical center manifold M^'^(r) at the R2 bifurcation can be parametrized locally by 
(r, as 

(68) u^uo+CiVi+^2V2 + H{T,0, re [0,2r], e = (a,6)eM', 
where H satisfies H{2T,^) — _ff(0,^) and has the Taylor expansion 

(69) H{t, = h^oei + hn^i^2 + ^/io2?|+ 

D 2 2 

where all functions hij are 2T-periodic, the eigenfunction corresponding to eigenvalue —1 is given by 

0, re [o,r], 

0, 
0, 



0, re [0,T], 

0, 

0, 



-V2{t) for T e [0,T]. 

The functions vi and V2 exist due to Lemma 5 of |23j . The functions hij of (69) can be found by solving 
appropriate BVPs, assuming that (2) restricted to VF'^(r) has normal form (9). As in the generalized period- 
doubling case, we first deduce a property for these functions hij. More general than in the GPD case, we 
here have that u(r, ^i, ^2) = u(t + T, —^1, — ^2)- This implies that 

and thus 

from which follows that hij{T + T) = hij{T) for i + j even and hij^r + T) = —hij^r) for i + j odd, for 
T S [0, T]. Taking these periodicity properties into account, we can reduce our observations to the discussion 
of the interval [0, T] instead of [0, 2T]. 

The coefficients a, a and b arise from the solvability conditions for the BVPs as integrals of scalar products 
over the interval [0, T]. Specifically, those scalar products involve among other things the quadratic and cubic 
terms of (3) near the periodic solution uq, the eigenfunction vi. The adjoint eigenfunction ip* associated to 
the trivial multiplier is the T-periodic solution of (29). The adjoint eigenfunction is the unique solution 
of the problem 

r vl{T)+A'^{r)vl = 0, re [0,r], 
(73) <^ vtiT) + vm = 

[ g{vlV2)dT-l = 0. 



{vi - A{t)vi = 

viiT)+vi{0) = 

Jo {vi,vi)dT -I = 

and the generalized eigenfunction by 

{V2 - A{t)v2 +Vi = 

V2{T)+V2{0) = 

Jo {v2,vi) dr = 

with 

(72) vi{t + T) := -vi{t) and V2iT + T) :-- 
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Note that we can indeed require this normahzation since V2 is the last generahzed eigenfunction of the 
original problem and therefore not orthogonal to all the eigenfunctions of the adjoint problem. We further 
define the generalized adjoint eigenfunction V2 as the unique solution of 

r i>2*(r) + AT(r)«2* ~vl = 0, T e [0,T], 

(74) <^ vUT)+vm - 

since, as above, is not orthogonal to V2- Moreover, we have 

^ ^^(^)) ^2, ^2)dr = 1^ {vt,V2)dr = 1. 

Note that 

(75) / {V2,Vi)dT= I {vl,Vi)dT= [ {V2,V2)dT = Q. 

Jo Jo Jo 

To derive the normal form coefficients, we proceed as in the previous sections, namely, we substitute (68) 
into (2), and use (3) as well as (9) and (69). 

By collecting the ^°-terms we get the identity 

Wo = F{uo). 

The linear terms provide two other identities, namely 

■ill = A{t)vi, V1+V2 = A{t)v2, 

in correspondancc with (70) and (71). 

Collecting the ^|-terms gives us an equation for ho2 

ho2 - A{T)ho2 = B{t; ^2,^2) - 2/iii, 

to be solved in the space of functions satisfying ho2{T) = /io2(0). In this space, the differential operator 
^ — A{t) is singular and its null-space is spanned by uo- The Fredholm solvability condition 



/ {ip*,BiT;V2,V2)-2hil)dT = 

Jo 



gives us a normalization condition for function hn, i-e. 



f 

Jo 



{<p* ,hii) dr =]- I {(p*,B{t;v2,V2)) dr. 



2 , 

By collecting the ^1^2-terms we obtain the differential equation for hn 

fin - A{T)hn = B{t; vi,V2) - /i20, 
which must be solved in the space of functions satisfying hniT) — /iii(O). The orthogonality condition 



T 

{(fi* ,B{t;vi,V2) - /120) rfr = 
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gives us a normalization condition for /120, i-e. 

(76) / {^*,h2o)dT= [ {<p*,B{t;vuV2)) dr. 

Jo Jo 

By collecting the ^i-terms we find an equation for /120 

(77) h2o - A{T)h2o = B{t; vi,v-i) - 2auo, 

to be solved in the space of functions satisfying h2o{T) = /i2o(0). In this space, the differential operator 
^ — A{t) is singular and its null-space is spanned by uq. The Predholm solvability condition 



T 

{tp* , B[t; vi,vi) — 2auo) dr = 



leads to the expression 



(78) {>p*,B{T;vi,Vi))dT, 

where vi is defined in (70). 

With a defined in this way we have to find a normalization condition which makes the solution of (77) 
unique. Indeed, if /120 is a solution of (77) with h2o{T) = /i2o(0)) also /120 = ^20 + £1^0 is a solution, since 
uo spans the kernel of the operator ^ — A{t) in the space of T-periodic functions. The projection along the 
space generated by iiQ is fixed by solvability condition (76). So /120 can be found as the unique solution of 
the EVP 

( h2o- A{T)h2o- B{T;vi,vi) + 2aF{uQ) = 0, Te[0,r], 

(79) I /i2o(r) -/i2o(0) = 0, 

I !o {^* Mo) dr = J^{(P*,B{t;vi,V2)) dr. 

In the line of the previous observations, we can define /in as the unique solution of the BVP 

( hii-A{T)hii-B{T;vuV2)+h2o = 0, Te[0,r] 

(80) I MT)-/in(0) = 0> ^ 

I !o ^'P* dT = IJ^{(p*,B{t;v2,V2)) dT, 

with /120 defined in (79). 

By collecting the ^^-terms we get an equation for /130 

(81) hso - A{T)h3o = C{t; vi,vi,vi) + 3S(t; wi, /120) - 6au2 - 6avi, 

which again must be solved in the space of functions satisfying /i3o(T') = — /i3o(0). Taking the integral 
condition of (73) into account, we obtain 

1 

(82) ^^6J +3-B(t;?;i,/i2o) - 6aA(r)i;i) dr, 

where a is defined by (78), /i2o is the solution of (79) and vi and vl are defined in (70) and (71), respectively. 
As remarked before, it is important to note that if /i3o is a solution of (81) with h^^iT) = /i3o(0), also 
hao = hso + CsqWi is a solution, since vi spans the nullspace of the operator ^ — A{t). 
Collecting the ^1^2-terms we get the equation for /121 

(83) A21 - A{T)h2i = -/130 - 26t;2 - 2av2 - 2avi + C{t; vi,vi,V2) + B{T;h2o,V2) + 2B{T;hii,vi), 
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to be solved in the space of functions satisfying h2i(T) = — /i2i(0). The solvabihty of this equation implies 
/ {vl,-h2,Q - 2bv2 - 2av2 - 'iavi + C{t; vi,vi,V2) + B{t; h2o,V2) + 2B{t; hn,vi))dT = 0. 

^0 

Notice that the 6^1^2-term in normal form (9j is resonant: in fact we cannot use the freedom on h^o to make 
the normal form parameter b zero since 

/ {vt,h3o)dT= / {v*,h3o + elQVi)dT = / {v*,h3o)dT. 
Jo Jo Jo 

Using the normalization from (73 ) and (75) gives us the following expression for b 

1 F 1 f"^ 

^=0 / {vi,-2aA{T)v2+C{T;vi,vi,V2)+B{T;h2o,V2) + '2B{T;hn,vi))dT- - / {vl,h3o)dT. 
^ Jo ^ Jo 

There is no need to compute explicitly the cubic expansion of the center manifold since the last term of this 
sum can be rewritten as 

+^^^^0 ''2*'^'°^'^^ = [^''*^' ("i +^^^0 ^''^'^'^ 

= ~ {v2,C{t;vi,vi,vi) + 3B{T;vi,h2o) — Gav2 — QaAvi)dT 
Jo 

= - {v2,C{t;vi,vi,vi) + 3B{T;vi,h2o) - QceAvi)dT, 
Jo 

obtaining 

1 

(84) b=- {vl,-2aA{T)v2 + C(r; vi,vi,V2) + B{t; /i2o, ^^2) + 2B(t; hn,vi))dT 

1 

+ / {v2,C{T;vi,vi,vi) + 3B{T;vi,h2o)-QaAvi)dT, 
^ Jo 

where /i2o is defined in (79) and a calculated in (78). Notice that, since h^o appears on the RHS of equation 
(83) , wc have two degrees of freedom on /121 • In fact, if /121 is a solution of (83), also /121 = /121 +£21^1 +^30^2 
is a solution since 

dh2i .7 dh2i , I f '^^a ^ \ dh2i , . / 

By collecting the ^i^|-tcrms we got the equation for hi2 

hi2 - A{T)hi2 = C{t,vi,V2,V2) + B{T,vi,ho2) + 2B{T,V2,hn) - 2/i2i, 

to be solved in the space of functions satisfying h\2{T) = — /ii2(0). The Fredholm solvability condition 
implies that 

/ {vI,C{t, Vi,V2, V2) + B{t, vi,ho2) + 2B{t, V2, hii) - 2h2i)dT = 0. 
Jo 

As mentioned before, /121 has a component in the direction of V2, which is not orthogonal to the adjoint 
eigenfunction vl, so it is possible to impose 

/ {vl,h2\)dT =1- I {vI,C{t,vi,V2,V2) + B{T,vi,ho2) + 2B{T,V2,hn))dT. 
Jo ^ Jo 
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This condition defines £|q uniquely; the freedom of £21 gives us as usual another freedom on hi2 in the 
direction of V2- 

Finally, collecting the ^f'terms gives 

ho3 - A(t)/io3 = C{t, V2,V2, V2) + 3B{v2, /102) - 3/112, 
to be solved in the space of functions satisfying ho3{T) ~ — /io3(0)- The Fredholm solvability condition is 

T 

{vI,C{t, V2,V2,V2) +iB{v2,ho2) ~ 3hi2)dT = 0, 



which can be satisfied imposing 



T 1 '■^ 







vl,hi2)dT ^ - / {vI,C{t,V2,V2,V2) +3B{v2,hQ2))dT. 
-J JQ 



This last condition determines the value of £21 ^nd thus the third order center manifold expansion is uniquely 
determined. However, since this third order expansion of the center manifold is not needed for the compu- 
tation of the critical coefficients, we don't write down those conditions. 

3.6 Strong resonance 1:3 bifurcation 

The three-dimensional critical center manifold VF^(r) at the R3 bifurcation can be parametrized locally by 
(r, ^) as 

(85) u = u„+^v + (v{r) + H{t, ^, f"), r G [0, 3T], ^ e C, 

where the real function H satisfies H{3T,^,^) = H{0,^,^) and has the Taylor expansion 

HiT,i,o - lh2oe+hii^^+ho2^^ + h,oe + lh2ie^ 

(86) + ^/ii2ee' + J/iosf + 0(|en, 

with hij{ST) — hij{0) and hij — hji so that ha is real. The eigenfunction v is defined as the unique solution 
of the BVP 

r v{t)-A{t)v = 0, T e [0,T], 

(87) I v{T)--e'^v{0) = 0, 

I J^{v,v)dT-l = 0, 

and extended on the interval [0, 37"] using the equivariance property of the normal form, i.e. 

v{t + T) := e'^v{T) and v{t + 2T) e'-^v{T) for r e [0,r]. 

The definition of the conjugate eigenfunction v follows immediately. These functions exist due to Lemma 2 

of ng. 

As usual the functions hij can be found by solving appropriate BVPs, assuming that (2) restricted to 
VF^(r) has the periodic normal form (10). Also here we can deduce a property for the functions hij. The 
definition of v[t) in [0, ST] states that u(r, f) = u{t + T, e'^^'^/^^, e'^'^/^f). Therefore, 

k,i k,i 

and thus 
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for T e [0,r]. This for example implies that hkk is T-periodic. These periodicity properties allow us to just 
concentrate on the interval [0, T] 

The adjoint cigcnfunction Lp* corresponding to the trivial multiplier is the unique T-periodic solution of 
BVP (29). The adjoint eigenfunction v* satisfies 

r v*iT) + A^{T)v* = 0, T e [o,T], 

(88) I v*{T)-e'^v*iO) = 0, 

[ J^{v*,v)dT-l = 0. 

Similarly, we obtain v*. 

We now write down the homological equation and compare term by term. The constant and linear terms 
give us as usual 

uo = F{uo), v-A{t)v = 0, V - A{t)v = 0. 
From the or ^^-terms we obtain the following equation (or its complex conjugate) 

h2o - A{T)h2o = B{t; v, v) - 2bv, 

to be solved in the space of functions satisfying h2o{T) = e'^/i2o(0). In this space the operator ^ — A{t) 
has a range space with codimension one which is orthogonal to v* . So one Fredholm solvability condition is 
involved, namely 

{v*,B{t;v,v) - 2bv)dT = 0, 



which makes it possible to obtain the value of the coefficient b. In fact, 

(89) b^^j\v*,B{r;v,v))dT. 
Using this value for b we can find /120 as the unique solution of the BVP 



(90) 



h2o-A{T)h2o-B{T;v,v) + 2bv = 0, t € [0,T], 
/i2o(T)-e*^/i2o(0) = 0, 
Jq {v*,h2o)dT = 0. 



By collecting the ^^-terms we obtain an equation for hn 

hii - A(r)/iii = B{t; v, v) - aiiiQ, 

to be solved in the space of functions satisfying hi\{T) ~ /iii(O). The Fredholm solvability condition with 
if* gives us the value of a\ 

(91) ai = / {^*,B{T;v,v))dT. 

Jo 

With ai defined in this way, let hn be the unique solution of the BVP 

{till — A{T)hii — B{t]v,v) + aiiio = 0, Te[0,r], 
hii{T) - hii{0) = 0, 
J^{^\hll)dT = 0. 

Finally, collecting the ^^^-terms gives an equation for /121 

/121 - A{T)h2i = C(t; V, V, v) + 2B(t; v, hn) + B{t; v, /i2o) - 2cv - 2bho2 - 2aiV, 
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to be solved in the space of the functions satisfying h2i{T) = e^^i /i2i(0). Therefore, there must hold that 

/ {v*, C{t; V, V, v) + 2B{t; v, hn) + B{t; v, /120) - 2cw - 26/io2 - 2aiv)dT = 0, 
^0 



such that parameter c of (10) is determined by 



1 

c=- / {v* ,C{t]v,v,v) + 2B{T;v,hii) + B{T;v,h2o) - 2aiAv)dT, 
^ Jo 

where ai and b are defined by (91) and (89), respectively, and v, hn and /120 are the unique solutions of the 
EVPs (87), (92) and (90). 

Since we have to check the solvability of all the equations up to the maximal order of the normal form, 
we also collect the ^"^-terms and obtain 

/130 - ^(r)/i3o = C(t; V, v, v) + 3B{t; v, h2o) - 6bhn - 6a2Uo, 

to be solved in the space of the functions satisfying hso{T) = /i3o(0). Therefore, there must hold that 

/ (</?*, C(r; V, V, v) + 3B{t; v, /120) — Qbhn — 6a2Uo)dT = 0, 
^0 

which determines the value of a2, namely 

r-T 



ce2 = {(fi* ,C{t;v,v,v) + 3B{T;v,h2o))dT. 
Jo 



Remark that as in the Chenciner case v is not uniquely determined. Indeed, when v is a solution of (87) 
and 7 G C with 7^7 = 1, then 'yv is also a solution. Then the adjoint function is given by jv*, and b and 
/120 are replaced by 7^& and 7^/1201 respectively. The normal form coefficient c stays the same. However, the 
normal form coefficient b is multiplied with 7^. This doesn't affect the bifurcation analysis since this normal 
form coefficient only has to be different from zero, and obviously 7 7^ 0. Moreover, the analysis around the 
bifurcation point is independent from the sign of b. 

3.7 Strong resonance 1:4 bifurcation 

The three-dimensional critical center manifold W^(r) at the R4 bifurcation can be parametrized locally by 
(r, ^) as 

(93) u = uo + ^v + ^v{t)+H{t,^,^), re [0,4T], ^eC, 

where the real function H satisfies ff(4T, ^,^) = H{0,£^,^) and has the Taylor expansion 

H{t,^,^) = lh2oe + hii^^+^ho2^ + lh3oe + lh2ie^ 

(94) + Ih,2^e + Iho3e + O{\^f), 

with hij (4T) = hij (0) and hij = hji so that hu is real, while v is defined by 

v-A{t)v = 0, re [0,T], 

(95) { v{T)-e'iv{0) = 0, 

J^{v,v)dT-l = 0, 
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extended on [0,4r] using the equivariance property of the normal form, i.e. 

v{t + T) e*^u(T) = iw(r), 

v{t + 2T) := e"w(r) = -w(t), 
v{t + 3T) := e'^v{T) ^ -iv{T), 

for T e [o,r]. 

The definition of the conjugate v follows from this. These functions exist due to Lemma 2 of [23] ■ As 
usual the functions hij can be found by solving appropriate BVPs, assuming that (2) restricted to T/F^(r) 
has the periodic normal form (11). Similar to the Rl:3 case, there holds that 

/ifci(T) = /i,,(T + T)(e-"/2)'^-(e^-/2)', 

for r e [o,r]. 

The adjoint eigenfunction ip* is defined by the T-periodic solution of (29) and v* satisfies 

r v*{t) + A^{t)v* = 0, T e [o,r], 

(96) <^ v*{T)-e'iv*{0) = 0, 

Similarly, we obtain v*. 

The constant and the linear terms give the identities 

uo = F{uo), v~~A{t)^0, v-A{t)v^O. 

From the or f ^-terms the following equation (or its complex conjugate) is obtained 

(97) h2o-A{T)h2Q^B{T;v,v). 

Notice that this equation is non-singular in the space of functions satisfying h2o{T) = — /i2o(0). So /120 is 
given as the unique solution of the BVP 



(98) 



h2o-A{T)h2o~B{T;v,v) = 0, Te[0,r], 

/^20(r)+/l20(0) = 0. 



By collecting the ^^-terms we obtain an equation for hn 

/ill - A{T)hu = B{t; V, v) - aiUg, 
to be solved in the space of functions satisfying hii{T) = /iii(O). The Fredholm solvability condition 

T 

{(fi*, B{t; V, v) — aiiiQidr = 

gives us the possibility to obtain the value of ai, namely given by (91). With this value of ai, hn is the 
unique solution of BVP (92). 

The ^f^-terms give an equation for hi2 

hi2 - A{T)hi2 = C(t; v, v, v) + B{t; u, /102) + 2B(r; u, hn) - 2cv - 2aiv, 

to be solved in the space of functions satisfying hi2{T) — —1/112(0). The Fredholm solvability condition leads 
us to the value of c, namely 

(99) i^*^^i^-^^^^^^) + B{r;v,ho2) + 2B{T;v,hn)-2aiA{T)v)dT, 



3 Computation of critical coefficients 



27 



where ai is defined in (91), and v, hn and /102 are the unique solutions of the EVPs (95), (92) and the 
complex conjugate of (98). Taking the complex conjugate gives us the critical coefficient c. 
Finally, by collecting the -terms we obtain an equation for /103 

ho3 - A{T)ho3 = C{t; v, v, v) + 3B{t; v, ho2) - 6dv, 

to be solved in the space of the functions satisfying ho3{T) = i/io3(0). The non-trivial Fredholm solvability 
condition 

i-T 

{v*,C{t; V, V, v) + 3-B(r; v, /102) — Qdv)dT = 



/ 

Jo 



10 

gives us the value of the critical coefficient of (11), namely 

(100) d=l [ {v*,C{t; V, V, v) + 3B{t; v, /io2))dr. 

6 Jo 

So we finally obtain the value of 

which makes it possible to understand which bifurcation scenario of the R4 resonance we have. 

Also in this case v is not uniquely determined, since for every 7 € C with 7^7 = 1, 71; is a solution. 
Then the adjoint eigenfunction is given by 7U*, and h2o is replaced by 7^/120- The normal form coefficient c 

stays the same, but instead of d wc get ^"^d. However, this again doesn't influence the bifurcation analysis 
since the study is determined by the above defined a for which we need only \d\. 

3.8 Fold-Flip bifurcation 

The three-dimensional critical center manifold M^'^(r) at the LPPD bifurcation can be parametrized locally 
by (r, as 

(101) U = Uo + ^lVi+^2V2 + H{T,0, TG [0,2T], ^ = (^l,^2)eM^ 
where H satisfies H{2T,^) = H{Q,S) and has the Taylor expansion 

(102) H{t, = \h20il + + \ho2^l+ 

\hzoei + ^/i2i^f6 + \hMl + \hoi^l + 0(||ef ), 

while the eigenfunctions V\ and V2 are given by 

( vi- A{t)vi- F{uo) = 0, re[0,T], 

(103) I «i(r)-«i(o) = 0, 

i /(f (z;i,F(uo))dT = 0, 

and 

r V2-A{t)v2 = 0, re [o,r], 

(104) I V2{T)+V2{Q) = 0, 

[ /(f {V2,V2)dT -1 = 

with 

(105) Vi{t + T) := vi{t) and V2{t + T) := -V2{t) for r € [0,T]. 
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The functions vi and V2 exist because of Lemma 2 and Lemma 5 of (23j . The functions hij can be found 
by solving appropriate EVPs, assuming that (2) restricted to ly^(r) has the periodic normal form (12). 
Moreover, similar as before, u(t, ■^i, ^^2) — u{t + T, fi, —^2) such that 

for T e [0,T]. As before, we will reduce all computations to the interval [0,T]. 

The coefficients of the normal form arise from the solvability conditions for the EVPs as integrals of scalar 
products over the interval [0, T]. Specifically, those scalar products involve among other things the quadratic 
and cubic terms of (3) near the periodic solution uqi the generalized eigenfunction vi and eigenfunction ^2, 
and the adjoint eigenfunctions Lp* , Vi and V2 as solution of the problems 

r ^*+A^[t)^* = 0, re [0,T], 
(106) I ^*{T)-^*{0) = 0, 

[ {ip*,V,)dT-l = 0, 



vi + A^{T)vi + ip* = 0, Te[o,r], 

(107) { vl{T)-vm = 0, 



{vl,vi)dT = 0, 



and 



v*2+A'^{t)vI = 0, re [o,r], 
(108) { v*^{T)+vm = 0, 

J^{vlV2)dT-l = 0. 

Note that the integral conditions are possible due to the spectral assumptions at the LPPD point. The 
following orthogonality conditions hold automatically 



(109) / {^\F{uo))dT= / {^*,V2)dT= / {vlv2)dT 

Jo Jo 

{v*2,V,)dT^ [ {v;,F{uo))dT -0, 



and since we have normalized the adjoint eigenfunction associated to multiplier 1 with the last generalized 
eigenfunction, we have for free that 



(110) 1 = 1^ {^*,v,)dr ^- A^(r)^ «i*,«i)dr 



{vl, ^— - A(T)j v,)dr = {vlFiuo))dT. 

As usual, to derive the normal form coefficients we write down the homological equation and compare 
term by term. 

By collecting the constant and linear terms we get the identities 

iiQ = F{ua), vi = A{t)vi + F{uo), V2 = A{t)v2, 

and the complex conjugate of the last equation. 

By collecting the ^^"terms we find an equation for /120 

(111) /120 - A{T)h2o = B{t; vi,vi) - 2a2oWi - 2a20'«o + 2wi, 
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to be solved in the space of functions satisfying h2o{T) = /i2o(0)- In this space, the differential operator 
^ — A{t) is singular and its null-space is spanned by uq. The Predholm solvability condition 



i 



T 

{ifi*, B{t; vi,vi) - 2a2oWi - 2q!2oWo + 2ui) dr = 



10 

gives us the possibility to calculate parameter 020 of our normal form, i.e. 

(112) {<fi*,B{T;vi,Vi) + 2A{T)vi) dr. 

With 020 tuned in this way equation (111) is solvable, for any value of parameter a2o- As in the Cusp 
of cycles case, we are free to choose parameter a2o as we want, and we take a20 = 0. This choice will not 
influence our final conclusion about the kind of situation we are in. 

In order to make the solution of (111) unique, we have to fix the projection on the null-space of the 
operator, more specifically in the direction of F{uq). Therefore we impose the orthogonality with the adjoint 
generalized eigenfunction , and obtain /120 as the unique solution of the BVP 

r h2o-A{T)h2o-B{T;vi,vi) + 2a2oVi+2a2oF{uo)-2A{T)vi-2F{uo) = 0, Te[0,T], 

(113) <^ h2o{T) - h2o{0) = 0, 

By collecting the ^1^2-terms we obtain a singular equation for hu 

/ill - A{T)hn = B{t; vi,V2) - 611^2 + V2, 
to be solved in the space of the functions that satisfy hu{T) = — /iii(O). The Predholm solvability condition 



T 

{v2 ,B{t; vi,V2) - biiV2 + V2) dr = 



10 

gives us the possibility, using (104) and (108), to calculate coefficient 611 of our normal form 

(114) 611=/ {v*2,B{t;vi,V2)+A{t)v2) dr. 

Jo 

With 611 defined in this way we can compute hn as the unique solution of the BVP 

( hu-A{T)hn-B{T;vi,V2) + hiV2-A{T)v2 = 0, T€[0,r], 

(115) I hu{T) + hu{0) = 0, 
[ J^{v*2,hn)dr = 0. 

Collecting the ^f-terms gives a singular equation for ho2 

(116) ho2 - A{T)ho2 = B{t; ^2,^2) - 2ao2f 1 - 2ao2Wo 

and, since this equation has to be solvable, the following Predholm solvability condition is involved 

/ {(p*,B{t;v2,V2) - 2ao2Vi - 2ao2Wo) dr = 0, 

JO 

from which we obtain an equation for ao2 

1 

(117) "'02=2y {'fi*,B{T;V2,V2)) dr. 
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3.8 Fold-Flip bifurcation 



So (116) is solvable, for any value of the parameter ao2- For simplicity, we take ao2 = 0. 

Notice that also here, the solution of (116) is orthogonal to the adjoint eigenfunction t^*. Since we have to 

fix the projection in the direction of eigenfunction mq, we define /io2 as the unique solution of 

( ho2- A{T)ho2- B{T;v2,V2)+2ao2Vi+2ao2F{uo) = 0, re[0,T], 

(118) I ho2{T) - ho2{0) = 0, 

[ J^{vt,ho2)dT = 0. 

By collecting the ^f-terms we get a singular equation for hso 

hso - A(r)/i3o = C(t; vi_,vi,vi) + SB{t; /120, ^^i) - 6020/120 - 6030^1 + 3/i20 - 6a3oWo - 6020*1, 
where the Fredholm solvability condition 

/ {ip*,C{T;vi,vi,vi) + 3B{T;h2o,vi) - 6020/120 - ^a^oVi + 3/i2o - 6a3oUo - 6a20^'i) dT = 
Jo 

gives us the value of 030 
1 

(119) aso = g / {^* , C{t;vi,vi,vi) + 3B{T;h2o,vi) -6a2oh2o 

+3{A{T)h2o + B{t;vi,vi)) + 6{1 - a2o)AiT)vi) dr- oao- 

Similarly, by collecting the ^1^2-terms we get a singular equation for /121 

h2i - A{T)h2i = C{t; vi,vi,V2) + B{t; ft20,W2) + 2B{T;hn,vi) - 2a2ohn 
- 26ii/iii - 2621^2 + 2/111 - 2Q2oi'2, 

where the Fredholm solvability condition 

/ (^2 , C{t;vi,vi,V2) + B{T;h2o,V2) + 2B{T;hii,vi) 
Jo 

-2a2o/iii - 2611/111 - 2621W2 + 2/iii - 20201)2) dr = 

produces the value of 621 

(120) 

1 

&21 = 2 / {v2,C{t;vi,vi,V2) + B{T;h2o,V2) + 2B{T;hii,vi) -2a2ohii 

-2611/111 + 2(A(r)/iii + B{t; vi,V2)) + 2(1 - a2o)A{T)v2) dr - bn. 

By collecting the ^i^|-terms we obtain a singular equation for hi2 

hi2 - A{T)hi2 = C{t;vi,V2,V2) + B{T;ho2,vi) + 2B{t; hii,V2) 
-26ii/io2 - 2ao2/i20 - 2ai2Wi + /102 - 2ai2Wo - 2ao2i'i, 

where its solvability requires the following equality 

/ {ip* , C{t;vi,V2,V2) + B{T;ho2,vi) +2B{T;hii,V2) 
Jo 

-2611/102 - 2ao2/i20 - 2ai2Wi + /102 - 2ai2Uo - 2Q;o2i'i) dr = 0, 
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such that 
(121) 



1 



ai2 = - / {(P*,C{t;vi,V2,V2) + B{T;hQ2,vi) + 2B{T;hii,V2) - 2biiho2 

-2ao2/i20 + A{T)ho2 + B{t; W2, ^2) - 2ao2^(T)ui) dr - 002- 
Finally, the ^|-ternis give us the value of the last needed critical coefficient. The equation for Hq^ is 
^03 - ^(t)/io3 = C{t; V2,V2,V2) + 3B{t; /102, "2) - 6002/111 - 6603^2 - 6q!o2W2, 
with Fredholm solvability condition 



(wj, C(r; U2, V2, V2) + 5B{t; /102, W2) - 6002^111 - 6603^2 - 6ao2W2) dr = 0, 



and thus 
(122) 



bo3 = 7: 



{v2,C{t;v2,V2,V2) + 3S(t; /io2,W2) -6ao2^ii - 6ao2v4(r)u2) dr. 



4 Implementation issues 

Numerical implementation of the formulas derived in the preceding sections requires the evaluation of in- 
tegrals of scalar functions over [0,T] and the solution of nonsingular linear BVPs with integral constraints. 
Such tasks can be carried out with continuation software such as auto |l4j, CONTENT [31], and matcont 
[nllTO]. In these software packages, periodic solutions to (1) are computed with the method of orthogonal 
collocation with piecewise polynomials applied to properly formulated BVPs. 

The standard EVP for the periodic solutions is formulated on the unit interval [0, 1], so that the period 
T becomes a parameter, and it involves an integral phase condition: 



(123) 



i(T)-T/(x(T),a) = 0, re [0,1], 
x{0)-x{l) = 0, 
Jo{^ir),iir)) dr = 0, 



where ^ is a previously calculated periodic solution to a nearby problem, rescaled to [0, 1]. 

In the orthogonal collocation method problem (123) is replaced by the following discretization: 



(124) 



3=0 \j=0 



Xo,Q — a;Ar_i,; 



= 0, 
= 0, 



Af-1 m-1 

'^i,j{^i,3^ii,3) + '^N,o{xn,0,£,N,q) = 0. 

1=0 j=Q 

The points Xij form the approximation of x{t) with m + 1 equidistant mesh points 

j 

Ti,j = H {n+i - Ti), j = 0, 1, . . . , m. 



in each of the N intervals [Ti,ri_|_i], where 



— To < Ti < ■ ■ ■ < tn — 1. 
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4.1 Discretization symbols 



The £ij (r)'s are the Lagrange basis polynomials, while the points (j — 1, . . . , m) are Gauss points 
[9], i.e. the roots of the Legendre polynomial of degree m, all relative to the interval [rj,Ti+i]. 

With this choice of collocation points Qj, the approximation error at the mesh points has order of 
accuracy to, 

where h = maxi=i,2,...,Ar{ii}, U — Ti — Ti_i (i = 1, . . . , N), while for the coarse mesh points the error has 
order of accuracy 2to, 

\\x{n) - x,,o\\ ^ Oih^"') 

( "superconvergence" ) . 

The integration weight of Tij is given by Wj+ii^+i for < i < N — 1 and < j < m. For 
I = 0, ... ,7V - 2, the integration weight of Tj_,„ (r^^m = n+i^o) is given by crj,,„ = Wm+iU+i + wiU+2, and 
the integration weights of tq and rjv are given by witi and Wm+itiy, respectively. In the above expressions, 
Wj+i is the Lagrange quadrature coefficient. 

The numerical continuation of solutions of (124) leads to structured, sparse linear systems, which in auto 
[13] and CONTENT 31j are solved by an efficient, specially adapted elimination algorithm that computes the 
multipliers as a by-product, without explicitly using the Poincare map. To detect codim 1 bifurcations, one 
can specify test functions that are based on computing multipliers |13[ 114] or on solving appropriate bordered 
linear EVPs [I5]. 

4.1 Discretization symbols 

It is convenient to discretize all computed functions using the same mesh as in (124). For a given vector 
function rj G C^([0, 1],R"') we consider three different discretizations: 

• rjM G ]ij(Wm+i)n^ ^Yie vector of the function values at the mesh points; 

• rjc £ K^™"^ the vector of the function values at the collocation points; 

• Vw — [ vwl ] G M^™" X M" , where rjwi is the vector of the function values at the collocation points 
multiplied by the Gauss-Legendre weights and the lengths of the corresponding mesh intervals, and 
r]w2 = '7(0). 

Formally we also introduce the structured sparse matrix Lcxm that converts a vector t/m of function 
values at the mesh points into a vector rjc of its values at the collocation points, namely, rjc = LcxmVm- 
This matrix is never formed explicitly; its entries are approximated by the £ij(Ci_fc)-coefficients in (124). 
We also need a matrix Acxm such that AcxmVm = {A-{t)'r]{t))c ■ Again this matrix need not be formed 
explicitly. On the other hand, we do need the matrix {D — TA{t))cxM explicitly; it is defined by {D — 
TA{t))cxMriM = {ilit) - TA(t)r]{t))c- Finally, let the tensors BcxMxm and CcxMxMxm be defined by 
BcxMxMriiMmM = iB{t;T]iit),T]2it)))c and 

CcxMxMxMViMmMmM = (C(i; ?7i (t) , 7^2 (i) , ??3 (i)))c 

for all r]i S C^([0, 1],E"). (These tensors are not formed explicitly.) 

Let f{t),g(t) € C''([0,1],M) be two scalar functions. Then the integral f(t)dt is represented by 
EiIo^E"li^j(/c)ij<j+i = Eilo^ where {fc)i,j = fiQ.j) and Wj is the Gauss-Legendre 

quadrature coefficient. The integral f{t)g{t)dt is approximated with Gauss-Legendre by fy/^gc ~ fwi-^CxMQM, 
where equality holds if g{t) is a piecewise polynomial of degree to or less on the given mesh. For vector 
functions f{t),g{t) G C°([0, 1], R"), the integral {f{t),g{t)) dt is formally approximated by the same ex- 
pression: fy/^9c ~ fwi^CxM9M, where again we have equality if g{t) is a piecewise polynomial of degree to 
or less on the given mesh. Concerning the accuracy of the quadrature formulas, we first note that accuracy 
is not an important issue for the phase integral in (123), as this equation only selects a specific solution 
from the continuum of solutions obtained by phase shifts. Similarly, the discretization of the normalization 
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integrals does not affect the inherent accuracy, including superconvergence at the main mesh points Tj of the 
solution of the discretized BVP. Discretization of integrals, as specified above, follows the standard Gauss 
quadrature error, which has order of accuracy 2m if, as mentioned, the function g{t) is a picccwisc polynomial 
of degree m or less on the given mesh and if f{t) is sufficiently smooth (in a piecewise sense). Otherwise, 
still assuming sufficient piecewise smoothness, the order of accuracy of the numerical integrals is m + 1 if m 
is odd, and to + 2 if m is even. In particular, for the often used choice m = 4, the integrals would then have 
order of accuracy 6. 



4.2 Cusp of cycles bifurcation 

The first task is to rescale the computed functions to the interval [0, 1]. We start by defining ui{t) = uoiTt) = 
uo{t) for t e [0, 1]. The linear EVP's (15) and (21) are replaced by 



(125) 

with v{t) = V\{t/T) and 



i>i(i) - TA{t)v^{t) - TF{ui{t)) = 0, i e [0, 1], 
wi(l)-vi(0) = 0, 
^l{v^{t),F{u^{t)))dt = 0, 



ipt{t)+TA^{t)ipm 



0, [0,1], 

0, 



J^{v^l{t)Mt))dt-i = 0, 



where <p*(r) = ip\{T/T)/T, respectively. We also need to rescale the adjoint generalized eigenfunction defined 
by (22): 

vl{t) + TA^{t)vl{t) + Tiflit) = 0, te [0, 1], 

vt{l)-vl{0) = 0, 
{vl{t),v,{t))dt = 0, 

with v*{t) = vI{t/T)/T. Now, ai = and h2 i is the unique solution of the BVP 
(126) 

/i2,i(0 - TA{t)h2,i{t) - TB{t; vi{t),vi{t)) - 2TA{t)vi{t) - 2TF{ui{t)) + 2aiTF{ui{t)) 

/i2,i(l)-/i2.i(0) 
!l{vl{t)M,i{t)) dt 

where h^ir) = h2,i{T/T). Therefore, we obtain 



0, t e [0,1], 

0, 
0, 



(127) 



^ , {^l{t),-6aiA{t)vi{t) +3A{t)h2,i{t) + 'iB{t:vi{t),vi{t)) 
6 Jo 

+6A{t)vi{t) + 3B{t;h2,i{t),vi{t)) + C{t;vi{t),vi{t),vi{t))) dt. 



We now determine the matrix solutions for the several functions. We compute v-^m by solving the 
discretization of (125) 



(128) 



{D - TA{t))cxM 
6q - di 



p 





Vim 
a 



Tfc 

Oraxl 





where a equals zero since the M x M upper left part of the big matrix is singular, g{t) = F{ui{t)), and p is 
obtained by solving the following system 



[P^ b] 



{D-TA{t))cxM 
So - Si 

r2 



ri 



= [0 



Mxl 
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4.3 Generalized period-doubling bifurcation 



with ri and r2 random vectors. In the solution of this system 6 = 0; in (128 ) we then use the normahzed p. 
This technique guarantees that we are working with well-defined systems. 

We will compute (Pi^y instead of (/J^m since ^ii^ can be computed with a matrix which is very similar to 
the matrix from (128). Formally, the computation of y^^y is based on Proposition C.l from the appendix 
since 

" ^l{t)+TA^{t)^l{t) 

^^(l)-^^(o) 

" D-TA{t) 
6o - 5i 



we have that 



^1 



is orthogonal to the range of 



{^1)1 



{D-TA{t))cxM 
So - Si 



= 0, 

. By discretization we obtain 
= 0. 



Therefore, (filyir can be obtained by solving 



w 



{D-TA{t))cxM 
So - Si 



P 



q' 

where a equals zero and q is the normalized right null- vector of 



= [0 



Mxl 



. We then approximate 



iD-TAit))cxM 
So - Si 

I = {ipl{t),vi{t))dt by Ii = {fD'^^LcxMViM- fiw is then rescaled to ensure that Ii = 1. 

It is more efficient to compute vly^r than w^ji^, since will be used only to compute integrals of the form 
/o < '^^i(*)>C(*) > dt. Proposition C.5 learns us how to determine v^. Indeed, 

ii-TAh 
h{0) - h{l) 

for all appropriate functions h, thus f * can be obtained by solving 







_ vm 





]y-{ 






■ h ' 













[{vD 



T 
W 



{D-TA{t))cxM Vie 

So — Si Onxl 




[{TifiDwi^CxM 



0]. 



where a equals zero and p is defined above. 

Next, (/i2,i)m is found by solving the discretization of (126), namely, 



{D-TAit))cxM 
So - Si 



P 




i^Dwi^CxM 

with a = and p defined above. 

Finally, (127) is approximated by 

1 ' .*^T 



h2,lM 

a 



TBcxMxmVimVim + '^TAcxmVim + 2T£(c 

Onxl 





C = -{'fi)wi{'^-^CxMh2,lM + ^BcxMxMVimVim 

+&AcxMVim + SBcxMxMh2,lMVlM + CcxMxMxMVimVimV\m)- 



4.3 Generalized period-doubling bifurcation 

As done in Section 4.2, we first rescale the computed quantities to the interval [0, 1]. The linear BVPs (28) 
and (29) are replaced by 

vi{t)-TA{t)vi{t) = 0, [0,1], 
vi{l) + vi{Q) = 0, 
Jl{vi{t),vi{t))dt-l = 0, 
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where v{t) = vi{t/T)/VT, and 

ri{t) + TA^{t)^l{t) = 0, [0,1], 

^1(1) -^1(0) = 0, 

J^{^t{t),F{u^mdt-l = 0, 

with ip*{t) = ip\{T/T)/T. This leads to the expression 

1 

ai.i = 2 {^l{t),B{t;v^{t),vi{t))) dt, 

with ai,i = Tai. 

The adjoint eigenfunction, defined by (30), determines the following rescaled v^: 

r vt{t) + TA^it)vi{t) = 0, [0,1], 

(129) <^ ^;i*(l)+^;i*(0) = 0, 

I Io{^t{t)Mt))dt-l = 0, 

with v*iT) =vI{t/T)/VT. 

Let /i2,i be the unique solution of the BVP 

/i2,i(t) - TA{t)h2,i{t) - TB{t;vi{t),vi{t)) + 2ai,iTF(wi(t)) = 0, t € [0, 1], 

h2,i{l) - h2,i{0) = 0, 
/o((^t(t),/i2,i(t))dt = 0, 

where /i2(t) = /i2,i(r/T)/r, and /i3,i be the unique solution of the BVP 

/i3,i(0-2'^(t)^3,i(t)-T(7(t;wi(t),Wi(i),t;i(t))-3TB(t;wi(t),/i2,i(i)) + 6ai,iTA(f)wi(t) = 0, t€ [0,1], 

hsAl) + h3,m = 0, 

/oKW,/j3,iW) = 0, 

where hsir) = /i3,i (r/T)/(VTT). 
Therefore, we obtain 

"2,1 = TTT / (<pt(t),£»(t;t;i(i),t;i(0,^^iW,^^i(0) + 6C(*;^^i(*)>^^i(0./i2,i(i)) + 3B(t;/i2,i(t),/i2,i(i))+ 

4B(t;t;i(i),/i3,i(t))-12ai,i(^(t)/i2,i(i) + S(t;wi(i).^^i (*)))) + af,i> 

with a2,i = r^Q!2. 

Now, /i4,i is obtained as unique solution of the following BVP 

/i4,i (i) - TA{t)hi,i{t) - TD{t; «i (i), «i (t) , V, (t) , (t)) - 6rC(i; t;i (t) , (t), /12.1 W) 
-3TS(t; /i2,i(i), /i2,i - 4TS(t; ?;i(t), /i3,i(i)) + 12ai,ir(^(t)/i2,i(i) + vi{t),Vi{t)) 

-2ai,iF(ui(i))) + 24a2,iTF(wi(t)) = 0, t e [0, 1], 
/i4,i(l)- Vi(0) = 0, 
J^{^t{t),h4,i{t))dt - 0, 

with /i4(r) = /i4a(T/T)/T2. 

Finally, we can write down the critical normal form coefficient 

1 

{vl {t) ,E{t;vi (i) , vi (t) , t;i (i) , vi (t) , t;i (t) ) + 10£l(i ; t;i (i) , t;i (t) , vi (i) , /i2,i (t) ) 



120 r2 

+ 15C(i; /i2,i /i2,i(i)) + 10C(t; /i3,i W) + 10B(t; /i2,i /i3,i W) 

+ 5B{t;vi{t),h4,i{t)) - 120a2,iA{t)vi{t) - 20ai,iA{t)h3^i{t)) dt. 
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4.3 Generalized period-doubling bifurcation 



We now come to the implementation details in MatCont. We compute vim by solving 



{D-TA{t))cxM 
So + Si 



T 



with pi and gi the rescaled unit vectors of 











Ocxi 


Pi 




Vim 




Oraxl 







ai 




1 



{D-TA{t))cxM 
6o + Si 



T 
^2 











Ocxi 


n 








Onxl 











1 



and 



[pj 0,3] 



{D-TA{t))cxM 
60 + Si 



n 




= [Om 



xl 



where ri and r2 are random vectors. Every is equal to zero. We normalize Vim by requiring X^^g^ Sj^o '^j((^iM)i,j, (^^im) 
1, where aj is the Lagrange quadrature coefficient. 

The discretization of (129) can be computed with the same matrix, see Proposition C.2 of the appendix, 



[{v*i)l a] 



{D-TA{t))cxM 
So + Si 



T 

3i 



Pi 




= [0 



Mxl 



where a = 0. We then approximate I = {vl{t),vi{t))dt by Ii = {vD'^^LcxmVim- 'f^iv^ then rescaled 
to ensure that /i = 1. 

An analogous matrix is used to compute ifl^y: 



w 



{D-TA{t))cxM 
So - Si 



P 




[Omxi 1] J 



where q is the normalized right null- vector and p the normalized left null-vector of 



{D-TAit))cxM 
So - Si 



and a is equal to zero. In what follows we will use these definitions for p, q,pi and qi. We then approximate 
/ = {ip\{t),F{ui(t)))dt by Ii = (}Pi)y^^^gc and normalize V'im/ ^'^ ensure that Ii = 1. 
This then leads to the discretization of the expression for ai^i: 



(130) 



Q!l,l = -{f*i)wiBcxMxMVlMVlM- 



Now, /i2,i and h^^i are found by solving the following systems: 

h 



{D-TA{t))cxM ^ 
So - Si P 



'-2,1M 

a 



TBcxMxmvimvim — ^(Xi,iTgc 

Onxl 




and 



{D-TA{t))cxM 

So + Si 

{vDwi^CxM 



Pi 





^3,1M 

b 



rhs 

Onxl 
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respectively, with 

rhs = TCcxMxMxMVimVimVim + ^TBcxMxMVlMh2,lM — Qoii^iTAcxMVlM- 

a and b will be zero. 
Thus, 



(131) a2,l = -^i^Dwi y^CxMxMxMxMVlMVlMVlMVlM + QCcxMxMxMV\MV\Mh2,lM 

+ 3-BcxMxM^2,lM^2,lM + ^BcxMxMV\Mhs,\M 

— 12ai,i(AcxM/l2,lM + BcxMxMVimVim) 



"1,1- 



Then, /14 i is found by 



{D-TA{t))cxM ^ 
So - Si 

i^Dwi^CxM 



a 



rhs 

Onxl 





with 



rhs = TDcxMxMxMxMVimVimVimVim + QTCcxMxMxMVlMV\Mh2,lM 
+3TBcxMxMh2,lMh2,lM + ^TBcxMxMVlMha^lM 

-12q!i,iT(Acxm/i2,im + BcxMxmVimVim - 2q!i,i5c) - 24a2,iTfifc 

and a = 0. 

Now, we have all ingredients for the computation of the normal form coefficient 

(132) e = ]^20T'2 ('^l)w"i ^-^^CxMxMxMxMxM^^lMflM^^lMl^lM^^lM + ^0DcxMxMxMxMVlMVlMVlMh2,lM 
+ 15CcxMxMxMVlMh2,lMh2,lM + ^OCcxMxMxMVlMVlAlhs^iM + 10-BcxMxM^2,lM^3,lM 

+ 5BcxMxMVlMh4,lM — 120a2,l AcxMt^lM — 20ai,i Acxm/13,1M • 



4.4 Chenciner bifurcation 

As before, we rescale the computed quantities to the interval [0, 1]. The linear BVPs (41), (29) and (42) are 
replaced by 



(133) 

with v{t) = Vi{t/T)/VT, 



vi (t) - TA{t)vi (t) + iojT vi{t) = 0, i e [0, 1] , 
z;i(l)-t.i(0) = 0, 
S^{vi{t),vi{t))dt-l = 0, 



^l{t)+TA^{t)^l{t) = 0, [0,1], 
j^{ipl{t),F{ui{t)))dt-l = 0, 
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4.4 Chenciner bifurcation 



with V3*(r) = ^pI{t/T)/T, and 

vl{t)+TA^{t)vl{t)+iujTvl{t) = 0, ie[0,l], 
vl{l)-vm = 0, 
^l{vl{t),v^{t))dt-l = 0, 

where v*{t) = vI{t/T)/\/T, respectively. 
Then /120 is approximated by 

(.-iA^ j h20,i{t)-TA{t)h20,i{t) + 2iujTh20,i{t)-TB{t;vi{t),vi(t)) = 0, t G [0,1], 

^ ^ I /i2o,i(l)- Vi(0) = 0, 

where h2o{T) = h2o,i{T/T)/T . 

Before being able to compute the approximation to hw we need ai, defined by (44): 

"1,1= / {vl{t),B{t;vr{t),Vt{t)))dt, 
Jo 

with ai,i = Tai. Therefore, we get 

hnAt)-TA{t)hnAt)-TB{t;vi{t),vi{t)) + ai.iTF{ui{t)) = 0, ie[0,l], 

hn,i{^) - hn,i{0) = 0, 
{ipl{t),hiiAt))dt = 0, 

with hn{T)=hnAr/T)/T. 

Now, we can compute normal form coefficient c: 

i 

ci = --y {v*^{t),C{t-,viit),vi{t),vi{t))+2B{t;vi{t),hnAt))+B{t\Mi),h^^^^ 
with ci = Tc. With ci defined in this way, /121M can be computed as follows 

h2iAt) - TA{t)h2i,i{t) + i^Th2iAt) - TC{t;vi{t),vi{t),vi{t)) - 2TB{t;vi{t),hn.i{t)) 

-TB{t-,h2oA't^),vi{t)) + 2iciTvi{t) + 2ai,iT{A{t)vi{t)-iuvi{t)) = 0, ie[0,l], 

/I2i,i(l)-/J2i,i(0) = 0, 

{vl{t),h2i,i{t))dt = 0, 

where /i2i(t) = h2i,i{r/T)/{VTT). 
Next, the rescaling of /iso gives us 

hsoAt) - TA{t)h3oAt) + ^^^Th3oAt) - TC{t-,vi{t),vi{t),vi{t)) - 3TB{t-,vi{t),h2o,i{t)) = 0, t £ [0,1], 

/I30,l(l)-/l30,l(0) = 0, 

with /i3o(t) = h3o,i{T/T)/{VTT). 

Now, we need the rescaled h^i 1 before being able to compute coefficient a2 1 
(135) 

ft3i.i(t) -T^(t)/?.3i.i(i) +2«wT/i3i.i(t) -TZ?(t;wi(t),wi(t).t.i(t),5i(t)) 
3rC(i; ^;i(t), wi(t), /iii,i(t)) - 3TC(t; «i(t), 5i(t), /^20,i W) - 5TB{t; /iii,i(t), /i20.i(i)) 

-3TB(t; z;i(i), /i2i,i W) ' ^^^(i; «i W, ^iso.iW) 
+6icir/i2o,i(t) + 3ai,iT{A{t)h2o,i{t) - 2iw/i2o,i(i) + B{t; vi{t), vi{t))) = 0, t G [0, 1], 

/i3i,i(l)-/i3i,i(0) = 0, 
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where /i3i(t) = h^i,i{T /T) /T^ , so 



a2,i = I / {^l{t),D{t- vi{t), viit), vi{t), vi{t)) + C{t; vi{t), vi{t), hm,i{t)) + 4C(t; vi{t),vi{t), /iiia(t)) 

+ C{t;vi{t),vi{t), ViW) + 2S(t; W) + 2S(t; /ii2,i W) + ^(i; /i02,i /i20,i W) 

+ 2B(i; /i2i,i (0) - 4ai,i(A(i)/in,i(i) + vi{t), vi{t))))dt + al^, 

with a2,i = T'^a2- 

Now, we still need h22,i{t)- 

{ h22,i{t) - TA{t)h22,i{t) - TD{t; v,(t).v,{t), vi{t),vi{t)) - TC{t; vi{t), vi{t),ho2.i{t)) 
4TC{t; h^iAt)) - TCit- v,{t), .^(t), ^20.1 W) - 2rB(t; /in.iW) 

^2TB{U vi{t),hi2.i{-t)) - TB{t: /io2.i /i20.i (*)) - 2TB{t; /i2i,i(i)) 
+4ai,ir(A(t)/in,i(i) + S(t;t;i(<),l;i(i))-ai,iF(ui(t))) + 4a2,i7'^("iW) = 0, te[0,l], 

/l22,l(l)-/l22,l(0) = 0, 
J^{ipt{t),h22,l{t))dt = 0, 

where /i22(r) = h22,i(.T/T)/T'^. 

Therefore, we can compute the critical coefficient e: 



12r2 



{vl (t) , (t) , (i), vi (t) , (t), vi (t)) + D{t; Vi (t) , Vi (t), (t) , /io2,i W) 

+6Z?(t; «i (t) , «i (i), wi (t) ,hns{t))+ ^D{t; vi (t) , vi (t) , wi (t), /i2o,i (0) 
+6C(t; /iiiaW, + 3C(t; /ii2,i(t)) 

+3C{t-,vi{t),ho2,i{t),h2oAt)) + 6C{t;vi{t),hii,i{t),h2o,i{t)) + 6C{t-,vi{^^ 
+C{t;vi{t),vi{t),h3oAt)) 

+3B{t; hi2,i{t), /i20,i(i)) + 6B(t; /i2i,i(t)) + 3B(t; /i22,i(t)) 

/1024 ^304 W) + 2S(t; «i (t), /1314 W) ^ 12a24^W«i W 
-6aisiAit)h2iAt) + 2B{t;viit),hii,iit)) + Cit-,viit),viit),viit)) 

+B{t;h2o,i{t),vi{t)) - 2ai^iA{t)vi{t)))dt + Q2,i«^ + ai,ii^ - al ^^i^. 

We now impose the matrix solutions for the functions. We compute Vim by solving the discretization of 
(133) 



{D-TA{t) + iu:TL)c,,M 
5o - 5i 



P2 





Vl,lM 

a 



Ocxi 

Onxl 
1 



with a = 0, where q2 is the normalized right null- vector of the complex matrix K 



{D-TA{t) + iuTL)cxM 
So - Si 



and P2 the normalized right null- vector of K^. This vector is then rescaled so that < Vi{t), Vi{t) > dt = 1. 
For the computation of (piM, we use Proposition C.l to obtain 



iD-TAit))cxM 
So - Si 



P 




= [Om> 



where a equals zero. We then approximate / = (ipl{t),F{ui{t))}dt by Ii = {^pDy^^LcxMgM and we 
rescale (fl^ so that /i = 1. 



40 



4.4 Chenciner bifurcation 



For the computation of vl wc apply Proposition C.3 from the Appendix. Since vl lies in the kernel of the 



vm 



over there defined (1)2, we have that 
by solving 

{D-TA{t) + iojTL)c>,M 



_L (/)i(C^([0, 1], C")). The eigenfunction vl is thus computed 



So — 5i 

12 



P2 




[0 



Mxl 



1 . 



We then approximate / = {v\{t),vi{t))dt by /i = (wJ)^^Lcxm^^im and we rescale v^y^ so that /i = 1. 
(134) is approximated by 



{D-TA{t) + 2ii^TL)c^M 
5q - 5i 



'20, IM 



TBcxMxMVim1>1M 



The coefHcient ai^i can be approximated as 

0^1, 1 = {^i)wi^CxMxMVimVim 

which gives then all the information to determine the real function hn^i: 



{D-TA{t))cxM 

So -Si P 

{'PDwi^CxM 



hii.iM 
a 



TBcxMxmVimVim — OL\^\Tgc 




where a equals zero. 

Then an approximation for the rescaled normal form coefficient ci is given by 

Cl = — -(t^i)vv-i(^^CxMxMxMWlM^^lMt'lM+2i?CxMxMt'lM/ill,lM+ScxMxMt^ 

Next, we determine the third order coefficients of the center manifold expansion, namely 



{D-TA{t) + iujTL)cxM 
(^i)iyi-^cxM 



P2 





/l21,lM 

a 



rhs 

Orixl 





where 



rhs = TCcxMxMxmvimvimvim + 2TBcxMxMViMhn,iM + TBcxMxMh20,iMViM — '^iciTLcxMViM 

-'^ai,lT{AcxMVlM - ii^LcxMVlM) 



and a = 0, and 



{D-TAit)+3iujTL)cxM 
So - Si 



30, IM = 



TCcxMxMxMVimVimVim + 3TBcxMxMVlMh20,lM 
Onxl 



The approximation to (135 ) is given by 



{D-TA{t) + 2iujTL)cxM 
So - Si 



IM 



rhs 

Onxl 



with 



rhs = TDcxMxMxMxMVimViMVimViM + 'iTCcxMxMxMVlMVimhu^iM + 3TCcxMxMxMVlMVlMh20,lM 
+3TBcxMxMhii^lMh20,lM + 3TBcxMxMVlMh2l,lM + TBcxMxMVimhzo^lM 

—6iciTLcxMh2o,iM — 3ai^iT{AcxMh2o,iM — 2ia;I/cxM/i20,iM + -BcxMxm^^im^^im) 
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while 



1 



Oi2,l = ^(Vl)wi i^CxMxMxMxMVlMVlMVlMVlM+CcxMxMxMVlMVlMho2,lM+^CcxMxMxMVlMVlMhll^lM 

+ CcxMxMxMVlMVlMh20,lM + 2ScxMxM/lll,lM/ill,lM + 2ScxMxM^^1m/112,1M + -BcxMxm/102,1m/120,1M 

+ '^BcxMxMVlMh21,lM — 4ai,i(^C'xM/ill,lM + BcxMxMVimVIm)) + i- 



Now, we are able to compute /i22,i by 



{D-TAit))cxM 
6o - Si 

ivDwi^CxM 



p 





h22,lM 

a 



rhs 

Onxl 





with 



rhs = TDcxMxMxMxmVimVimVimVim + TCcxMxMxMViMViMho2,iM — '^a2,iTgc 

+iTCcxMxMxMVlMVlMhii^lM + 2Ti?CxMxM/lll,lM/lll,lM + 2Ti?CxMxMWlM^12,lM 
+TCcxMxMxMVlMVlMh20,lM + TBcxMxMho2,lMh20,lM + '^TBcxMxMVlMh21,lM 

— 4ai,iT(AcxM/iii,iM + BcxMxmVimVim — aigc) 
and a = 0. For the fifth order coefficient of the normal form, we then obtain 



e = 



Y^ij^{''^l)wii^CxMxMxMxMxMVlMVlMVlMVlMVlM + DcxMxMxMxMV\MVlMVlMhQ2,lM 

+6DcxMxMxMxMVlMVlMVlMhii^lM + 'iDcxMxMxMxMV\MVlMV\Mh20,lM 
+6CcxMxMxMVim/ii1,1m/ii1,1M + 5CcxMxMxMVlMVlMhi2SM 

+3CcxMxMxMVlMho2,lMh20,lM + 6CcxMxMxM^^lM/lll,lM/i20,lM + QCcxMxMxMVlMVlMh21,lM 
+CcxMxMxMVlMVlMh30,lM + 3-BcxMxM/ll2,lMft20,lM 
+6BcxMxMhn,lMh21,lM + 3BcxMxMVlMh22,lM 

+BcxMxMho2,lMh30,lM + 2-BcxMxM?^1m/131,1M — 12a2,l^CxMWlM 

— 6ai^l{AcxMh21,lM + ^BcxMxMVluhii^iM + CcxMxMxMVimVimVim 



. Ci 



+BcxMxMh2Q,iMViM " 2ai,i Acxm^^im)) + a2,i«^ + ai^iij^ - 
4.5 Strong resonance 1:1 bifurcation 

Again, we rescale the computed quantities to the interval [0, 1]. The linear BVPs (56) and (57) are replaced 
by 

vi,i{t) - TA{t)vr,x{t) - TF{ui{t)) = 0, t £ [0, 1], 
t^i,i(l)-t^i,i(0) = 0, 
{vi,i{t),F{ui{t)))dt = 0, 



where vi{t) = wi^i (r/T), and 



V2,i{t) - TA{t)v2,i{t) + Tvi,iit) = 0, t€ [0, 1], 
t^2,i(l)-t^2,i(0) = 0, 
{v2,i{t),F{Mt)))dt = 0, 



with V2{t) = V2,i{t/T), respectively. 
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4.5 Strong resonance 1:1 bifurcation 



The rescaling of the adjoint cigcnfunction, defined by (60), gives us 

^l{t) + TA^{t)^l{t) = 0, ie [0,1], 

= 0, 

Jo M (0,^^2,1 = 0, 

where <^*(t) = ip\{t/T)/T . This leads to the expression 



2A{t)vi^i{t) + B{t- vi^i{t)))dt. 



Definition (61) of the first generalized adjoint eigenfunction is rescaled as 

vl^it) + TA^{t)vl^{t) - Tiplit) = 0, i e [0, 1], 
vlii^)-vlM = 0, 
/o Ki(*)>«2,i(i))rft = 0, 

with u*(t) = vI i{t/T)/T, which then gives all the information we need to compute the critical coefficient b 

(136) b= [ {ipl{t),B{t;vi,i{t),V2,i{t)) + A{t)v2,i{t))dt+ [ {vli{t),2A{t)vi,i{t) + B{t;vi,i{t),vi,i{t)))dt. 
Jo Jo 



The generalized eigenfunctions are computed as 

{D-TA{t))cxM „ 



and 



So - Si 

9Wi^CxM 

[D - TA{t))cxM 
So - ^1 

9Wi^CxM 



a 



P 




V2,1M 

a 



Tgc 

Onxl 




-Tvi,ic 

Onxl 





with p the left null-vector, as before. 

As always, we compute the adjoint eigenfunctions by means of the transposed of the usual matrix, and 
obtain here (filyi^ instead of 'filJ^f■ Formally, the computation of ipl^r is based on Proposition C.l from the 
appendix and as in the cusp of cycles case wc obtain then 



w 



{D-TA{t))cxM 
So - Si 



p 





[Omxi 1] 



where a equals zero. We approximate / = Jq {ipl{t),V2,iit))dt by Ii = (<Pi)f4^j-£'CxM^^2,iM and rescale ip^^^ 
to ensure that Ii ^ 1. 

Having found vi^im and pl^/^ the first normal form coefficient of interest can be computed as 

1 

2^ 



^ = 7;{y^l)wii^^CxMVl,lM + BcxMxMVi^imVi,im)- 



Now we still need i- The technique of Proposition C.5 from the appendix is used to obtain (wjf i)w, 
namely 

" {D-TA{t))cxM V2,IC " 

(5o — <5i 0„xi 
c?T 



[Ki) 



w 



[-T{^*l)lr^LcxM 0], 



where a = 0. 

Finally (136) is approximated by 

^ = {fDwi {BcxMxMVi^IMV2,1M + AcxMV2,1m) + (2^CxM^^l,lM + BcxMxMVi,IMVi,Im)- 
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4.6 Strong resonance 1:2 bifurcation 

As before, we rescale the computed quantities to the interval [0,1]. The linear EVPs (70) and (71) are 
replaced by 

r vi,i{t)-TA{t)vi,i{t) = 0, te [0,1], 
(137) <^ fi,i(l) + ?;i,i(0) = 0, 

with vi{t) = vi^i{T/T)/Vf, and 

V2,i{t) - TA{t)v2At) + Tvi^iit) = 0, t e [0, 1], 
f2,i(l) + «2,i(0) = 0, 

{V2,l{t),Vi,i{t))dt = 0, 

where V2{t) = V2,i{t/T)/Vt. 

The rescaling of the adjoint eigenfunction gives us 

iplit) + TA^(t)^l{t) = 0, [0,1], 
^l{l)-ipl{0) = 0, 

j^{^tit),Fiu,mdt-i = 0, 

where 'P*{t) = ip\{T /T) /T , so we can compute 

1 



with ai = Ta. With ai defined in this way, let /120 1 be the unique solution of the EVP 
(138) 

h2o,iit)-TA{t)h2o,iit)-TB{t;vi,i{t),vi,i{t)) + 2aiTF{ui{t)) = 0, te[0,l], 

h2o,i{l) - h2o,i{0) = 0, 
/o((^t(t),/i2o,i(t)) = J^{<fl{t),B{vi,i{t),V2,im dt, 

where /i2o('r) = h2o,i{T/T)/T. 

The rescaling of the adjoint eigenfunction corresponding with multiplier —1 and the adjoint generalized 
eigenfunction gives us 

vliit) + TA'^it)vl,it) = 0, t e [0,1], 

IoKiit)^^2At))dt-l = 0, 
with vI{t) = vI-^{t/T)/Vt, and 



v*2At) + TA^{t)v*2At) - Tvl,{t) = 0,t€ [0, 1], 

«2*l(l)+^'2*l(0) = 0, 

IoMt),viiit)) dt = 0, 



where V2{t) — W2.i(T/T)/\/r, respectively. 
The rescaled critical coefficient is then 



1 

(139) ai = - (v^i(t),C(t;wi,i(t),wi,i(i),t;i,i(t))+3B(i;t;i,i(t),/i2o,i(t))-6aiA(t)wi,i(i)) dt, 
6 Jo 

with ai = Ta. 
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4.6 Strong resonance 1:2 bifurcation 



We replace (80) by 



hn,i{t) - TA{t)hn,i{i) ' TB{t;vi,i{t), V2,i{t)) + Th2,y,{t) 

/iii,i(l)-/iii,i(0) 
!l{^l{t)Mi,i{t)) dt 



0, t€ [0,1], 

0, 

i^'(<^^(^),S(^;2,l(^),^;2,l(^)))rf^, 



with hii{T) = hii^i{T/T)/T, to finally obtain that 
1 

b=— {vl^it),-2aiA{t)v2,iit) + C{t;vi,i{t),vi,i{t),V2,i{t)) 
1 

+ {v2,i{t),C{t;vi,i{t),vi,i{t),vi,i{t))+3B{t;viAt),h2oAt)) - GaiA{t)vi,i 



{t)}dt. 



We now come to the second part of the implementation details and define the matrix solutions for the 
several functions. We compute wi,im by solving the discretization of (137) 



(140) 



{D-TA{t))cxM 
So + Si 

T 



Pi 











Ocxl 




Wl,lM 




Onxl 
1 




a 





where a = 0. We then normalize Wi,im by requiring J^i^^o^ Sjlo (i^i,iM)i,j) = 1, where aj is 

the Lagrange quadrature coefficient. 

We have now obtained the value of vi^i in the mesh points. However, since is used in the integral 
condition for V2,i, we have to transfer this vector to the collocation points and multiply it with the Gauss- 
Legendre weights and the lenghts of the corresponding intervals, giving us vector {vi^Wi- Then, V2,im can 
be found by solving the following system 



Pi 



{D-TA{t))cxM 

So + Si 
{vi,i)wj_LcxM 



V2,1M 

a 



-Tvi^ic 

Onxl 




where a equals zero. 

The adjoint eigenfunction corresponding to the trivial eigenvalue is computed with an analogous matrix 
as in (140): 



{D-TA{t))cxM 
So - Si 



p 





[0 



Mxl 



1 



while the adjoint eigenfunction corresponding to eigenvalue —1 can be found by 

{D-TA{t))cxM 



So + Si 



T 

Qi 



Pi 




[Om> 



1 



where ai and 02 are equal to zero. <plyy is then rescaled to make sure that {'Pi)'^^Lcxm9m = 1 and vl 

so that {vIA'^^LcxmV2,im = 1- 

Making use of Proposition C.6, we obtain the adjoint generalized eigenfunction by solving 



[{V2,i)w a] 



{D -TA{t))cxM V2,1C 
So + Sl Onxl 

ql 



[-T{vli)^^LcxM 0]. 
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Having found flyir^ and Vi^im, oli can be computed as 



Now, /i20,iM is found by solving the discretization of (138 ), namely, 



(D-TA(i))cxM ^ 
and (139) is approximated by 



h20,lM 

a 



TBcxMxmVi^imvi^im — "^aiTgc 

Onxl 

{V*)'wi^CxMxMVi^1MV2,1M 



dl = g(^^l,l)vFi {CcxMxMxMVi^imVi^imVi^im + 3BcxMxMt^l,lM/i20,lM — 6ai Acxmi^i.im) 



Next, hii,iM is found by 



(£) - TA{t))cxM 
So - Si 



i'P*)wiLcx 



M 



P 





^11, IM 

a 



TBcxMxMVi,imV2,1M — Th20,lC 
Onxl 

2{'Pi)'wi^CxMxMV2,1MV2,1M 



Finally, we obtain 

b= -^{vl^i)wi{-'^0liAcxMV2,lM 



+ CcxMxMxMVi,1MVi,1MV2,1M + -BcxMxm/120,1M1'2,1M + ^BcxMxMhn,lMVl,lM) 

+ ^(^2,1)^1 {(^CxMxMxMVi^imVi^imVi^im + 3-BcxMxM1'1,1m/120,1M — QOiiAcxMVl,lM)- 

4.7 Strong resonance 1:3 bifurcation 

As before, wc rcscalc the computed quantities to the interval [0, 1]. The BVPs for the eigenfunction and its 
adjoint belonging to eigenvalue are replaced by 



(141) 

with v{t) = Vi {t/T)/Vt, and 



vi{t)-TA{t)vi{t) = 0, t e [0,1], 
wi(l) - e*^wi(0) = 0, 
{vi{t),vi{t))dt-l = 0, 



vt{t) + TA'^itX^it) = 0, te [0, 1], 

vm-e'^-^vtiO) = 0, 
J^{vl{t),viit))dt-l = 0, 

where v*{t) = vI{t/T)/\/T. The rescaling of the adjoint eigenfunction corresponding to the trivial multiplier 
gives 



(142) 

with (p*(t) = ipliT/T)/T. 



^lit)+TA'^{t)ipl{t) = 0, ie [0,1], 



j^{^t{t),F{u,mdt- 



0, 
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4.7 Strong resonance 1:3 bifurcation 



These eigenfunctions make it already possible to compute the following 2 rescaled normal form coefficients 
(143) 

where an = Tai, and 



"1,1= / {ipl{t),B{vi{t),vi{t)))dt, 
Jo 



with bi = VTb. 

The rescaled second order functions in the center manifold expansion are solutions of 



h20,i{t) - TA{t)h2o,i{t) - TB{vi{t), vi{t)) + 2hTvi{t) 

/J20,i(l)-e''^/i20,i(0) 



with /i2o(r) = /i20,i(t/T)/T, and 
(144) 



hii,i{t)-TA{t)hii,i-TB{vi{t),Vi{t))+ai,iTF{ui{t)) 

/iii,i(l)-/iii,i(0) 



0, t€[0,T], 

0, 

0, 



0, [0,1], 

0, 

0, 



with /iii(t) = (r/T)/T. This all results then in 



c = 



2T 



{vl{t),C{vi{t), vi{t),vi{t)) + 2B{vi{t),hn,iit)) + B{vi{t), /i20,i(t)) - 2ai,iAvi{t))dt. 



We now come to the implementation details in MatCont. Eigenfunction vi, determined by (141), is 
computed by 

{D-TA{t))cxM 
So — e~*"f (5i 

with a = 0. Wc then normalize vim by requiring '^f^o^ '^J^q crj{{viM)i,j, {''^iM)i,j) = Ij where aj is the 

" {D-TA{t))cxM ' 



P3 




Vim 
a 




Ocxl 
Onxl 









1 



Lagrange quadrature coefficient. (73 is the normalized right null-vector of JsT = 
the normalized right null- vector of K^, with 9 = 

a 

Since wjf e Kei{^2), this function can be obtained by solving 

{D - TA{t))cxM 



So - e-^'Si 



and P3 



To compute the adjoint eigenfunction v^, we apply Proposition C.4 from the appendix with = 



(145) 



Sq — e ' 3 5i 
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[Omxi 1] 



v*y^r is rescaled such that {vD^^Lcxm^im = 1- The adjoint eigenfunction corresponding to eigenvalue 1 is 
discretized by 



(146) 



iD-TA(t))cxM 
So - Si 

T 



p 




[Omxi 1] 



where a equals zero, cplyy is rescaled so that {fDy^r^LcxAigM = 1- 



4 Implementation issues 



47 



The normal form coefficients ai^i and bi become 

(147) Q!i,l = {(fiDwiBcxMxMVlMVlM 

and 



bl = -{vD^^BcxMxMVimVim- 



By computing first the complex conjugate of we can use the same matrix as in (145), except for the 
last line which represents the integral condition, to get 



P3 



{D-TA{t))cxM 
(^r)m-^CxM 



^20, IM 

a 



TBcxMxmVimVim — 2biTvic 

Onxl 





To obtain the discretization of hu,i, the following system is solved 
{D-TA{t))cxM ^ 



(148) 



which gives us 



So - 5i 



i^Dwi^CxM 



hii^iM 
a 



TBcxMxmVimVim — ai,\Tgc 

Onxl 




C = ■^{vI)^^{CcxMxMxMViMVimViM + '2,BcxMxMVlMh\\^\M + BcxMxMVlMh20,lM — 2ai,l^CxMl'lM)- 

4.8 Strong resonance 1:4 bifurcation 

The eigenfunction and the adjoint eigenfunction corresponding to multiplier e*5 are given by the solution of 

vi(t)-TA{t)vi{t) = 0, t e [0,1], 
i;i(l) - e*tvi(0) = 0, 
/o (t;i(i),t;i(t))dt-l = 0, 

with v{t) = Vi {t/T)/Vt and 

vl{t)+TA'^{t)vl{t) = 0, te[0,l], 

vl{l) - e^^vliO) = 0, 

J^{vt{t),vi{t))dt-1 = 0, 

wherew*(T) = vI{t/T)/\/t, respectively. ip*,hii and ai are replaced by y^, /in, 1 andai,i, defined by (142), 
(144) and (143), respectively. 
The rescaling of (98) gives 

h2o,i{i)-TA{t)h2o,i{i)-TB{vi{t),vi{t)) = 0, te[0,l], 

/l20,l(l) + /l20,l(0) = 0, 

with /i2o(t) = /l20,l(T/T)/r. 

The critical coefficients are then given by 

1 



and 



6T 



{v*^ (i), C{vi {t),vi {t), vi {t)) + 3B{vi {t), /io2,i {t)))dt. 
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4.9 Fold-Flip bifurcation 



The eigenfunction and its adjoint, corresponding to the complex eigenvalue, are discretized by 



{D - TA{t))cxM 
So — e~'5(5i 



if 











Ocxl 


P3 




Vim 



















a 




1 



and 



{D - TA{t))cxM 
6o — e~'5 (5i 
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P3 





[0 



Mxl 



respectively. For the computation of the adjoint function we have applied Proposition C.4 with ^ = §. We 
then normalize vim by requiring X^il^o^ Sj=o '^j{{viM)i,j, {viM)i,j) = 1, where aj is the Lagrange quadrature 
coefficient, vl^/ is rcscalcd so that {vD^^Lc x mVim = 1- 

(146) and the corresponding normalization, (148) and (147) determine <^i,/iii,i and Q!i,i, respectively. 
An approximation to /120 is obtained by 



{D-TA{t))c>,M 
(5o + <5i 



'-20, IM 



TBcxMxMVimVim 
Onxl 



We are now able to compute the two needed normal form coefficients: 



1 



C = —{vl)yir^{CcxMxMxMVlMVlMVlM + BcxMxMVlMho2,lM + 2ScxMxMWim/111,1M — "^C^i^iAcxMVIm) 



and 



d = 1 (CcxMxMxMViMViMViM + 3i3cxMxMf 1m/i02,1m)- 



4.9 Fold-Flip bifurcation 

The rescaling of the eigenfunctions (103) and (104) gives us 



(149) 



and 



(150) 



i;i,i(t) - TA{t)vi,i{t) - TF{ui{t)) = 0, t G [0, 1], 
vi,i(l)-t'i,i(0) = 0, 
{vi,i{t),F{ui{t)))dt = 0, 



V2,i{t)-TAit)v2,i{t) = 0, t e [0,1], 

V2,l(l)+^'2,l(0) = 0, 

/o KiW, 1^2,1 W)cit-1 = 0, 



respectively, with vi{t) = vi^iir/T) and W2(t) = W2,i(t/T)/-\/T- 
The rescaled adjoint eigenfunctions can be obtained by solving 

r^{t)+TA^{t)^l{t) = 0, t€ [0,1], 
= 0, 

J^{vl{t),vi,i{t))dt-1 = 0, 



vl,{t) + TA'^{t)vl,it) + T^lit) = 0,tG [0, 1], 
<i(l)-<i(0) = 0, 
lo Kiit),viA{t))dt = 0, 
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v*,{t)+TA^{t)vl,{t) = 0, [0,1], 

^^2,l(l)+^^2,l(0) = 0, 

with ^*(t) = ipl{T/T)/T,vl{T) = vI^{t/T)/T and v*2{t) = vI^{t/T)/VT. 
The two coefficients in front of the ^^-terms are given by 



1 

^20 = 2 J {^*i{t),B{vi,i{t),vi,i{t)) + 2A{t)vi,i{t)) dt 



and a2o = 0. 

The second order functions of the center manifolfd expansion are defined by the EVPs 

f h2o,i{t)-TA{t)h2o,i{t)-TB{vi,i{t),vi,i{t))+2a2oTvi,i{t)+2a2QTF{uiit)) 

-2TA{t)vi,i{t) - 2TF{ui{t)) = 0, t e [0, 1], 

Vi(l) -^^204(0) = 0, 

lo Kiit)''^2o,i{t))dt = 0, 

with /120(t) = /l20,l(T/r), 

huAt) -TA{t)hnAt) -TB{vi,i{t),V2,iit)) + TbnV2,i{t) -TA{t)v2,iit) = 0, t € [0,1], 

ftii4(l) + 'iii,i(0) = 0, 
lo Ki(*),hii,i{t))dt = 0, 

with /iii(t) = /iii,i(t/T)/\/T and 

ho2At)-TA{t)ho2,i{t)-TB{v2,i{t),y2,i{t)) + '^ao2,iTvi,i{t)+2ao2,iTF{ui{t)) = 0, te[0,l], 

Vi(i)-Vi(o) = 0, 

/o Ki(*)./i02,i(t)>rft = 0, 
with /io2(t) = /io2,i (t/T)/T, where 

^11= / {v;i{t),B{vi,i{t),V2,l{t)) + A{t)V2,l{t)) dt, 

Jo 

1 

ao2,i = 2 (<^tW,BKi(i).«2,i(i))) rfi 

with ao2,i = Tao2 and ao2 = 0. 

The rescahng of the last four normal form coefficients of interest gives 

1 

030 = 7/ i'Pii't) : C{vi,i{t),vis{t),vis{t)) + 3B{h2o,i,vis{t)) - 6a2ah2o,i{t) 

+3{A{t)h2o,i{t) + B{viAt),vi,i{t))) + 6(1 - a2o)A{t)viAt)) dt - a2o, 



QJo 



1 

^21 = 7 / Ki(i).C'(^^i,i(*)>^^i,i(*),^^2,i(*)) + -B(/i20,iW,t^2,i(t)) + 2B(/in,i(t),vi,i(i))-2a2o/iii,i(i) 
Jo 

-26n/in,i(t) + 2(A(t)/in,i(t) + BiviAt),V2,i(t))) + 2(1 - a2o)A(t)i;2.i(t)) dt - 
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4.9 Fold-Flip bifurcation 



1 

ai2 = ^ ((pt(i),C(^^i,i(i).^^2,i(*),^^2,i(t)) +S(/io2,i(i)>^^i,i(*)) + 2S(/iii,i(t),t;2,i(i)) - 26ii/io2,i(i) 



-2ao2,i/i20,i(t) + A{t)ho2,i{t) + B{v2,i{t),V2,i{t)) - 2ao2,iA{t)vi,i{t)) dt - 



(102,1 



and 



1 

bo3 = ^ / {v2i{t),C{v2,i{t),V2,iit),V2,i{t))+3B{ho2,iit),V2,i{t)) - 6ao2,ihii,i{t) - 6aQ2,iA{t)v2,iit)) dt. 
The discretization of the functions determined in (149j and (150) is given by 



iD^TA{t))c.M „ 
60-61 ^ 

QWi^CxM 



ai 



Tgc 

Onxl 





and 



iD-TA{t))cxM 
So + Si 



T 

Qi 



Pi 




V2,1M 

a2 



Ocxl 
Onxl 
1 



with ai = 02 = 0. We normalize t;2iM by requiring X^^q^ Sjlo ^j((^2,iM)ij, (w2,iM)i,j) = 1, where Uj is 
the Lagrange quadrature coefficient. The implementation of the adjoint eigenfunctions is done by 



{D-TA{t))cxM 
60 - 61 

T 

q 



p 




[0 



Mxl 



{D-TA{t))cxM ivi,i)c 
^0 — f^l 0„xl 
flT 



= [T(<^f)T^LcxM 0] 



and 



[{V2,l)w «3] 



{D-TA{t))cxM 
60 + 61 



T 

Qi 



Pi 




= [Om> 



with oi = a2 = 03 = 0. <^i(4/ and V2^mr are then rescaled to ensure that (<Pi)^jicxMi'i,iM 

(^2,i)w'i-^CxMV2,1M = 1- 

The first needed normal form coefficient is given by 

^20 = 2^fl)wi {BcxMxMVi,1MVi,im + 2AcxMVi,im) ■ 

The second order functions can be obtained by solving 



1 and 



{D-TA{t))cxM 
60 - 61 
(^i,i)m-^c'xM 



P 




''20, IM 

a 



rhs 

Onxl 
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with 



rhs = TBcxaixmVi,imVi,im - '2a2oTvi^ic + ^TAcxmVi^im + 2Tgc, 
{D-TA{t))cxM 



Pi 



5n + h 



a 



rhs 

Onxl 





with 
and 



rhs = TBcxMxmVismV2sm - T6iiW2ac + TAcxmV2,im^ 
{D - TA{t))cxM 



P 



So - Si 

(^r,l)w'i'^CxM 



ho2,lM 

a 



rhs 

Onxl 





with 



rhs — TBcxMxmV2,imV2,im — 2ao2,iTwi jc"- 
The in these functions needed coefficients are given by 

bll — (l'2,l)w'i (-^CxA/xA/t^l,lMW2,lA/ + ^CxA/W2,lA/) 

and 



a02,l = 2(¥'l)vKi^C'xAfxA/W2,lA/W2,lA//- 



At last, we obtain 



130 — QifDwi {Ccx M xM X MViAMViAMViAM + 3-BcxAfx Af/l20,lA/'yi,lA/ — 6020^120, IC 
+ 3(A(7xAf/l20,lAf + -BcxAfxAfWl,lAfWl,lAf) + QAcxMVi^im) — O.20, 



&21 — 2 '-^2^1^W'i i^CxMxMxMVl,lMVlAMV2AM + -Bcx AfxA/^20,lMW2,lA/ + 2i?CxAfx A/'ill,lAfWl,lAf — 2a20^11,lC 

— 2&ii/lii^lC + 2(yl(7xM^ll,lAf + BcxMxMVi.imV2.1m) + '^AcxMV2,1m) — ^lli 



ai2 = ■^{fDwi {CcxMxMxMVi^imV2,1MV2,1M + BcxMxMho2,lMVlAM + 2BcxMxMhii^iMV2,lM — 26iift.02,lC 

— 2ao2,l/l20,lC + ^CxAf^02,lA/ + BcxMxMV2,1MV2,1m) 

and 

^03 = ^(^2,1)^^1 {CcxMxMxMV2,1MV2,1MV2,1M + 3i?Cx AfxA/^024A/^'2,lAf — 60024/11140) • 

5 Examples 

The computations in this section are done with matcont TIT. In particular, the bordering methods from 
|15| are used to continue the codim 1 bifurcations of limit cycles in two parameters. The algorithms described 
above for computing the normal form coefficients are also implemented in the current version of matcont. 
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5.1 Periodic predator-prey model 



5.1 Periodic predator-prey model 

Our first model is a periodically forced predator-prey model, studied in |34| , and described by the following 
differential equations 



(151) 



ep{x,t)y~dy, 



where x and y are the numbers of individuals respectively of prey and predator populations or suitable (but 
equivalent) measures of density or biomass. The parameters present in system (151) are the intrinsic growth 
rate r, the carrying capacity K, the efficiency e and the death rate d of the predator. The function p(x, t) is a 
functional response, for which the HoUing type II is choosen, with constant attack rate a and half saturation 
b{t) that varies periodically with as period one year, i.e. 

QjX 

p{x,t) = ^^^^ 6(t) = 60(1 -f e cos 27ri). 

Notice that this system can be made autonomous by adding the following two differential equations 

a = fj.a-ujp - {a^ + f3^)a, l3 ^ uja + fi^ - {a^ + 13^)13. 

In fact, if ^ is positive then this system has as asymptotic behavior a stable circular limit cycle of radius fi 
with angular velocity uj. Therefore, we can set fJ- — I and cj = 27r in order to obtain the forcing function 
sin 27rt as the flow of one of the variables of this subsystem with a particular phase shift that depends on the 
initial conditions. The system becomes 



(152) 



P = 27ra-H/3- (a2 +/32)/3. 



We have chosen this system as first example since it allows us to check if the computation of the ai 
normal form coefficients is correct. In fact, in a periodically forced system the time of the normal form 
should not depend on the coordinate, i.e. dr/dt = 1, and so all the ai coefficients must vanish. With fixed 
r ~ 2tt, K = e= l,a — Att and d = 27r we perform a bifurcation analysis in the remaining parameters 
(e, &o) obtaining the bifurcation diagram reported in Figure 1 Since the system is periodically forced, 
no equilibria are present. The blue curve, with label LPC2, is a limit point of cycles bifurcation curve of 
the second iterate, the magenta curves are supercritical Neimark-Sacker bifurcations (of the first or of the 
second iterate, respectively labeled with NSl and NS2) while the brown and green curves are period-doubling 
bifurcations, brown when subcritical and green when supercritical (with notation PDl, PD2, PD4 and PD8). 



We now analyze in detail all the detected codimension two points, reporting the scalar computed coeffi- 
cients explained in Section 4, 



5.1.1 The two GPD points 

In Figure 1 the LPC2 curve is tangent to the PDl curve in two different GPD points. In the first one, in 
(e, 60) = (0.319,0.412), the limit point of cycles curve is tangent to the subcritical period-doubling curve 
(type presented in Figure 12-(b)), while in the second one, in (e, 69) = (1-09, 0.218), the LPC2 curve is tangent 
to the supercritical part of the PD bifurcation curve (i.e. the type presented in Figure 12 -(a)). 
Performing the computation of the GPD normal form coefficients at the first point we obtain: 

• for the first equation of normal form (6) the two coefficients ai and a2, up to a scaling term T and 
computed through the formula (130) and (131) are zero, up to the accuracy of the computation. 
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02 0.4 0.6 OB 



Figure 1: Bifurcation diagram of limit cycles in model (152). Blue are limit point of cycles bifurcations, green 
period doubling bifurcations and magenta Neimark-Sacker birucations. Continue/dotted curves correspond 
to supercritical/subcritical bifurcations. 

• the normal form coefficient of the second equation, computed through formula (132) equals e = 
-58.2867. 

Notice that these results are in agreement with what we expected, i.e. that since the system is periodically 
forced the time doesn't depend on the distance from the critical limit cycle, and since we are in the case 
presented in Figure 12 -(b) the normal form coefficient e is negative. 

From the computation of the GPD normal form coefficients at the second critical point we obtain: 

• for the first equation of normal form (6) the two coefficients equal zero. 

• the normal form coefficient of the second equation has value e — 41.5442. 
Also in this case the obtained results are in agreement with the theory. 

5.1.2 The 1:1 and 1:2 resonance points 

We divide the 1:1 and 1:2 resonance points present in this model into two groups, namely the R2 point at 
(e, bo) = (0.337, 0.34) and the cascade of resonance points in the lower part of the graph. 
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5.1 Periodic predator-prey model 



The isolated R2 point forms the intersection of the NSl curve, the supercritical Neimark-S acker curve of a 
limit cycle with period approximately equal to 1, and PDl. The situation is thus the one depicted in Figure 
15 -(a). Performing the normal form coefficient computation we obtain: 

• for the first equation of normal form (9) holds that a — 0. 

• for the last equation of the normal form (9) we have (a, b) = (3.401426, —12.90745). 

Notice that the obtained results are in accordance with the theory (no secondary Neimark-Sacker curve 
implies that a > and supercritical Neimark-Sacker curve implies that & < 0). 

In the lower part of the bifurcation diagram a resonance cascade is present, which accumulates on the 
sequence of period-doubling curves. A zoom of this part is shown in Figure 2, Each resonance point of this 
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e 

Figure 2: The resonance cascade in model (152). In blue are limit point of cycles bifurcations, green 
period doubling bifurcations and magenta Neimark-Sacker birucations. Continue/dotted curves correspond 
to supercritical/subcritical bifurcations. 

cascade (except for the first Rl point on the LPC2 curve) can be seen in two ways: the first (and more natural 
way) is to see them as R2 points of the type represented in Figure 15 -(b) (so with a < and the sign of 
b dependent on the criticality of the incoming Neimark-Sacker curve), the second way is to see them as Rl 
points represented in Figure 14 -(b). Notice that the criticality of each NS curve of the cascade changes at 
the R2 point (as depicted in Figure 15 -(b)). 

As first general result we see that for the first equation of normal form (9) there holds that a = for 
all points (as expected since the system is periodically forced). We remark that for the computation of the 
normal form coefficients of the second equation, the tolerances have to be strong enough. The results are 
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The Rl point is in (e, 60) = (0.6044021,0.1769608). The period 2 hmit cycle Neimark-S acker curve 
NS2 starts tangentially to the LPC2 curve and is subcritical. In this situation we therefore expect the 
product of the two normal form coefficients of the last equation of (8) to be positive; the computed 
coefficients are (a, b) = (2.005489e - 6, 6.4806e + 9), such that ab = 1.299679e + 4. 

The R2 point is in (e,6o) = (0.7437713,0.1835935). On the left side of the R2 point the PD2 curve is 
supercritical, on the right side it is subcritical. The NS2 curve incoming in the R2 point is subcritical, 
while the NS4 curve outgoing at the R2 point is supercritical. Wc arc thus in the case depicted in 
Figure 15-(b) with time reversed. So we expect 6 > (subcritical incoming Neimarck-Sacker curve) 
and o < (there is an outgoing secondary Neimarck-Sacker curve). The computed coefficients at the 
R2 point arc {a.,b) = (—65.76676,16.26708). Notice that this point is also a degenerate Rl point for 
the NS4 curve. In fact, when we compute the normal form coefficient at this 1:1 resonance bifurcation 
point, we obtain (a, b) = (-1.113774823354237e-4, 6.116846870167980e+12). In this case the product 
ab = — 6.812790e+ 8 is negative, in accordance with the fact that the NS4 curve (for which we have an 
Rl bifurcation) is supercritical. 

The R2 point is in (e,6o) = (0.7427991,0.1861098). On the left side of the R2 point the PD4 curve is 
supercritical, on the right side it is subcritical. The NS4 curve incoming in the R2 point is supercritical, 
while the NS8 curve outgoing at the R2 point is subcritical. We are therefore in the case depicted 
in Figure 15 -(b). We expect 6 < (supercritical incoming Neimarck-Sacker curve) and a < (there 
is an outgoing secondary Neimarck-Sacker curve). The computed coefficients at the R2 point are 
(a, 6) = (-269.3681,-18.15061). 

The R2 point is in (e, 60) = (0.7439079, 0.1864190). On the left side of the R2 point the PD8 curve is su- 
percritical, on the right side it is subcritical. The NS8 curve incoming in the R2 point is subcritical: we 
are thus in the case depicted in Figure 15 -(b) with time reversed. Thus, we expect 6 > (subcritical in- 
coming Neimarck-Sacker curve) and a < (there is an outgoing secondary Neimarck-Sacker curve, since 
the cascade continues). The computed coefficients of the R2 point are {a,b) = (—921.7011, 16.58059). 

All the obtained results are in agreement with the theory. 

5.1.3 The 1:3 resonance points 

There are two 1:3 resonance points, one on NS2, the other one on NS4, as can be seen in Figure 2, These 
two points behave in a different way. The Neimark-S acker curve corresponding with the first point at 
(e, 60) = (0.709,0.179) is subcritical, so we expect 3?(c) to be positive. The Neimark-Sacker curve of the 
second point at (e, &o) = (0.743,0.185) is supercritical, so 3?(c) should be negative. To check whether we 
arc in a non degenerate case, we also have to look at b, however, as mentioned before, the sign of b is not 
relevant. We obtain 

• for the first R3 point we have that (6, 3?(c)) = (4.5567 - 4.4567i, 9.155003). 

• for the second R3 point we have that (6, 5i(c)) = (0.4049 + 12.1425i, -8.819864). 
These results are in accordance with the theory. 

5.1.4 The 1:4 resonance points 

There are two 1:4 resonance points, one on NS2, the other one on NS4, as can be seen in Figure 2, Also 
these two points behave in the same way as the 1:3 resonance bifurcation points. The Neimark-Sacker curve 
corresponding with the first point at (e, 60) = (0.675,0.177) is subcritical, so here we expect ^{A) to be 
positive. The Neimark-Sacker curve of the second point at (e,feo) = (0.743,0.185) is supercritical, so 3?(A) 
should be negative. Moreover, since those points are part of a resonance cascade, we should not have limit 
point bifurcations of non trivial equilibria, so we are in region I of Figure 17 . In order to assure that we are 
not in a degenerate case, we also need to check that d^O. We obtain 
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• for the first R4 point we have that (c, d) = (11.624-84.897i, 65 + 92.254i), so A ^ 0. 102999 -0.752278i. 

• for the second R4 point we have that (c, d) = (-8.5796-414.71i, -416.64-489. 17i), so A = -0.013352- 
0.645406i. 

So for both bifurcation points the value of A belongs to region I, and thus the results are in accordance with 
the theory. 

5.2 The Steinmetz-Larter model 

The following model of the peroxidase-oxidase reaction was studied by Steinmetz and Larter |41) : 

-kiABX - k3ABY + k7~ k-rA, 
-kiABX - k-iABY + ks, 
kiABX - 2k2X^ + 2k3ABY - k^X + ke, 
-k^ABY + 2k2X^ - k^Y, 

where A,B,X,Y are state variables and fci, k2, k^, k^, k^, fcg, fcy, fcg, and /c_7 are parameters. We fix all 
parameters as reported in the following table 



Par. 


Value 


Par. 


Value 


Par. 


Value 


Par. 


Value 


fci 


0.1631021 




1250 




0.046875 


fc4 


20 




1.104 


ke 


0.001 


fc-7 


0.1175 







and we perform a bifurcation analysis in the parameter space {kj, k^). A few curves are reported in Figure 
3, 

5.2.1 The 1:1 resonance points 

The two 1:1 resonance points have different nature, since in one the Neimark-Sacker curve rooted at the 
bifurcation point is supercritical, while in the other one it is subcritical. 
We obtain: 

• for the Rl point in (fc/jfeg) = (1.179554,0.7239571), the two coefficients of the last equation of (8) are 
equal to (a, 5) = (-0.003654362200739,0.735048055230916). Their product ab = -2.686132e - 3 is 
negative, which corresponds with the fact that the NS curve rooted at the Rl point is supercritical. 

• for the Rl point in {kT,ks) = (1.857676,0.9304220), the two coefhcients of the last equation of (8) 
are equal to (a, 6) = (-0.066429738171756,-2.156596806473489). Their product ab = 0.1432622 is 
positive, and indeed the NS curve rooted at the Rl point is subcritical. 

So we can conclude that the results are in accordance with the theory. 

5.2.2 The Chenciner points 

As can be seen in Figure 3 we have detected a CH point at (fcyjfcg) = (1.757356,0.9125773). The normal 
form coefficient at that bifurcation point equals 3?(e) — 1.391931, hence positive so the unfolding is the one 
depicted in Figure 13 -(b). In order to verify if the normal form computation is correct, one should use tori 
continuation techniques |37j, either recurring to Poincare maps 26, 8 or to the so-called invariance equation 
[m 1311 [321 [35]. However, these techniques are not stable, especially in critical cases like the one we have. 
In order to validate our result we thus have to do simulations. 

The obtained result is shown in Figure 4 , The indicated regions correspond with the regions from Figure 
13 The green curve between regions 2 and 3 is the supercritical Neimark-Sacker curve, the red one between 
region 1 and 2 is the subcritical Neimark-Sacker curve. For each point of the grid, we have started time 
integration from a point close to the orginal limit cycle (a 1 % perturbation) until an attractor was found. 



(153) 
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Figure 3: Bifurcation diagram of a limit cycle in model (153). In blue are limit point of cycles bifurca- 
tions, green period doubling bifurcations and magenta Neimark-Sacker birucations. Continue/dotted curves 
correspond to supercritical/subcritical bifurcations. 



The 1-norm of the X-coordinate of an orbit with time length 1000 along the attractor is shown in the 
colormap. In region 2 this attractor is the original limit cycle, in region 3 it is the inner torus arisen through 
the supercritical Neimark-Sacker curve. In region 1 the original limit cycle is unstable, and so the trajectory 
which starts nearby converges to another attractor. Between region 1 and 3 and region 1 and 2 happens 
a catastrophic bifurcation, i.e. a drastic change of the attractor, identified from the change of color which 
varies from blue to red. Right above the Chenciner point, the catastrophic bifurcation is the subcritical NS 
curve, while left below it is the limit point of tori (Tc) curve. Figure 4 shows that we obtain the scenario 
which corresponds with a positive second Lyapunov coefficient. 

5.3 The Lorenzl984 system 

This model, taken from [3^, is a meteorological model written by Lorenz in 1984 in order to describe the 
atmosphere. The equations of the model are 

{X = —y^ — — ax + aF, 
y ^ xy — bxz — y + G, 
z = bxy + xz — z, 

where {a,b,F,G) are parameters, with a = 0.25, & = 4. This model, as depicted in f40|, has most of the 
analyzed codimension two bifurcations of limit cycles. We report in Figure 5 a bifurcation diagram obtained 
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5.3 The Lorbnz1984 system 
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Figure 4: Simulations on a parameter grid (black points) of system (153). The magenta continue/dotted 
line is the supercritical/subbritical Neimark-Sacker curve. The color represents the value of the maximum 
of the first coordinate of the attractor reached through simulation from a point close to the limit cycle. A 
sketch of the state portrait is reported in Figure 13 -(b). 



with MATCONT in which bifurcations of equilibria (LP stands for limit point, H for Andronov-Hopf) are thicker 
and limit cycle bifurcations are thin. In particular, the blue curve is a limit point of cycles (LPC) bifurcation 
curve, the green ones are period-doubling (PD) bifurcations curves and the magenta ones are Neimark-Sacker 
(NS) curves. The codimension two points are marked with a red dot, and, as can be seen in the figure, almost 
all cases (except the Chenciner bifurcation and the fold-flip bifurcation) discussed in Section 3 are present in 
this model. In the sequel of this section we will investigate the normal form coefficients of each bifurcation. 

5.3.1 The Swallowtail bifurcation 

The first degeneracy we want to analyze is the vanishing of the coefficient c in the cusp of cycles normal 
form (5). This bifurcation, named Swallowtail bifurcation, is characterized, in our case, by the collision and 
disappearance of two cusp points of limit cycle. In order to get this codimension three bifurcation we analyze 
part of the blue dashed curve in Figure 5 for different parameter values of b. The result is shown in Figure 
6 : in the graph part of the limit point of cycles manifold is plotted in the {F, G')-plane for different values of 
parameter b G [2.91, 2.95] (from blue to red). In the table we can see the behavior of the critical normal form 
coefficient c, where it exists (the colors of the lines correspond with the bifurcation diagram) . Notice how 
the behavior of this codim 3 bifurcation is captured by a smooth vanishing of the normal form coefficient. 

5.3.2 The degenerate generalized period doubling bifurcation 

On the green curve PD4 of Figure 5 there arc two generalized period-doubling (GPD) points. In the first one 
the flip bifurcation curve passes from subcritical to supercritical, in the second one the opposite happens. 
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Figure 5: Bifurcation diagram of model (154). The thicker curves are bifurcation curves of equihbria, the 
thin curves are bifurcation curves of hmit cycles and invariant tori (blue, limit point of cycles, green period 
doubling and magenta Neiamrk-Sacker). Continue/dotted curves correspond to supercritical/subcritical 
bifurcations. 



Computing the normal form coefficient in the first case gives e = — 1.317656e — 3 < 0, therefore there is a 
limit point of cycles bifurcation curve that starts rightward tangent to the supercritical part of the period- 
doubling manifold, while in the second case e = 2.895460e — 3 > 0, and so the limit point of cycles bifurcation 
starts leftward tangent to the subcritical part of the PD curve. These conclusions can be seen in Figure 7, 
where the period-doubling curve is black, dotted when supercritical; in the upper panels are sketched the 
Poincare map of the limit cycle involved in the bifurcation. On the yellow curve the limit cycles sketched in 
green and red collide and disappear, while on the violet curve the two involved limit cycles are sketched in 
red and blue. 
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5.3 The Lorbnz1984 system 



5.3.3 The 1:1 resonance points 

Two Rl points are located on tlie LPC curve. Tliose two points should have different product of the normal 
form coefficients. In fact, in the first one (in {F, G) = (10.72, 0.522)) the Neimark-Sacker curve rooted at the 
bifurcation point is supercritical (i.e. the system has the behavior depicted in 14-(a)), while in the second 
one (in {F, G) = (9.81, 2.22)) the NS curve is subcritical (behavior similar to Figure 14-(b)). If we apply our 
analysis we obtain: 

• for the first Rl point (a, b) = (2.577, —1.2659), so the product ab = —3.26237 is negative. 

• for the second Rl point (a, b) = (—9.887, —2.005), so the product ab = 19.81858e is positive. 

These results are in accordance with the theory. The blue curve on Figure 8 is the limit point of cycles curve. 
The violet curves are the Neimark-Sacker curves of first iterate, the green curve is the Neimark-Sacker curve 
of second iterate. 

5.3.4 The 1:2 resonance points 

At the unique R2 point at {F,G) = (10.72,0.522) shown in Figure 5, the incoming Neimark-Sacker curve, 
namely NS, is subcritical (and so > 0), while the outgoing curve (which exists and so a < 0), namely NS2, is 
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Figure 6: Different limit point of cycles bifurcation curves in the {F, G)-plane for different values of the third 
parameter b. The parameter values are reported in the table. 
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Figure 7: Two generalized period-doubling points with different normal form coefficients on the period- 
doubling bifurcation curve PD4 of Figure 5 , 

supercritical, i.e. we are in the case reported in Figure 15 -(b) with time reversed. The coefficients computed 
at the 1:2 resonance point are (a, 6) — (—0.6329965,0.1785283), in accordance with the theory. 

At the Rl point located at (F, G) — (9.81, 2.22) starts a resonance cascade, as shown in Figure 8-(b). On 
that cascade we can find many resonance points which we will analyze in what follows. In particular, since 
the R2 points belong to a cascade they are of the same type as presented in Figure 15-(b) (so a < 0), with 
at each step a change of criticality of the incoming NS curve. As mentioned before, the first NS curve, born 
at the Rl point, is subcritical, so for the first R2 point we expect that 5 > 0, while for the second one 6 < 0. 
The obtained numerical results are 

• for the first R2 point at (F, G) = (9.9158,2.2978) we have that (a, fe) = (-1.3157,0.11076). 

• for the second R2 point at {F,G) = (9.9197,2.2978) we have that (a, 6) = (-2.6228,-0.0564). 
Results are in accordance with the theory. 
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5.3.5 The 1:3 resonance points 

We have several 1:3 resonance points at which we can have a closer look. There is one R3 point located 
on the NS curve and two R3 points on the NS2 curve. The R3 point corresponding to the first iterate 
is at {F,G) = (4.628,1.624), with a positive normal form cocfScicnt of the Nciniark-Sacker curve such 
that we are in the case represented in Figure 16 -(b). The R3 points corresponding to the second iterate 
are at {F,G) = (7.072,1.235) and {F,G) = (8.989,0.7394), where the Neimark-Sacker curve is both times 
supercritical, so we are in the situation depicted in Figure 16 -(a). 

• for the R3 point at (F, G) = (4.628, 1.624) we have that (6, SR(c)) = (0.1913 - 0.5464i, 6.185900e - 2) 

• for the R3 point at (F, G) = (7.072, 1.235) wc have that (6,3fJ(c)) = (-0.4461 -0.1901i, -3. 612302e-2). 

• for the R3 point at (F, G) = (8.989, 0.7394) we have that (6, JR(c)) = (-0.1285+0.0168i, -1.950822e-2). 

These results arc in accordance with the theory. 

There are also R3 points on the cascade (see Figure 8 -(b)). The first one corresponds with a subcritical 
NS curve, while the second one corresponds with a supercritical NS curve. 

• for the first R3 point, in (F, G) (9.8888, 2.2798) wc have that (6, 3fi(c)) = (-2.9582 - 0.3599i, 0.7383). 

• for the second R3 point, in (F, G) = (9.9187, 2.2978) we have that (6, K(c)) = (2.7447+3.539H, -0.3847). 
Also in this case all results are in accordance with the theory. 

5.3.6 The 1:4 resonance points 

There are 5 1:4 resonance points at which we will we have a look. One is located on the NS curve, two others 
on the NS2 curve and the last two lie on the resonance cascade (sec Figure 8 -(b)). We obtain 

• for the R4 point at (F, G) = (3.376, 1.647) we have that {c,d) = (0.0501 - 0.0746«, -0.0577- 0.5422i) 
and so ^ = 0.0918 - 0.1368i (subcritical NS curve, case I). 

• for the R4 point at (F, G) = (9.777, 0.595) we have that (c, d) = (-0.0151 - 0.1348i, -0.0266 - 0.0411i) 
and so ^ = -0.307829 - 2.752998i (supercritical NS curve, case VIII). 

• for the R4 point at (F, G) = (6.620, 1.390) we have that (c, d) = (-0.0417 - 0.9915i, -0.4283 - 1.0826i) 
and so A = —0.035841 — 0.851653i (supercritical NS curve, case I). 

For the first and the last point no further bifurcation analysis is possible to confirm the correctness of the 
results (present curves rooted at the point are global bifurcations of limit cycles). Instead it is possible 
to continue all local bifurcations of limit cycles rooted at the second R4 point, obtaining the result shown 
in Figure 9 Note that we haven't made the distinction between region VII and VIII since we have not 
computed the fold of torus curve typical for region VIII. 

For the cascade (see Figure 8-(b)) we have that the first point is on a subcritical NS curve, while the 
second one is on a supercritical NS curve. Moreover, since they lie on a cascade, they should be of type I. 

• for the first R4 point at (F, G) = (9.9159,2.2978) we have that (c, d) = (0.0518- 1.7633«, -2.0143 + 
0.4546j) and so A = 0.025102 - 0.853919i (subcritical NS curve, case I). 

• for the first R4 point at (F, G) = (9.9197, 2.2978) we have that (c, d) = (-0.0282 - 6.8152i, -10.8446 + 
2.1458j) and so A = -0.002550 - 0.616491i (supercritical NS curve, case I). 

Also in this case the results are in accordance with the theory. 



64 



5.4 The Extended Lorenz1984 system 




0.594 0.595 0.596 0.597 0.598 0.599 0.6 0.601 

G 



Figure 9: Bifurcation diagram at the R4 point at {F,G) — (9.777,0.595). In blue are the hmit point 
of cycles bifurcation curves, in violet the Neimark-Sacker curves. Continue/dotted curves correspond to 
supercritical/subcritical curves. 



5.4 The Extended Lorenzl984 system 

As done in i.33], it is possible to extend the Lorenzl984 system (154) by adding a fourth variable which takes 
the influence on the jet stream and the baroclinic waves of external parameters like the temperature of the 
sea surface into account. The obtained system is 



(155) 



X — —y — z — ax + aF 
y = xy ~ bxz ~ y + G, 
z = hxy + xz ~ z, 
u = —Su + jux + K. 



We use the parametervalues mentioned in [33], i.e. 

a = 0.25, 6 = 1, G = 0.2, 6 = 1.04, 7 = 0.987, F = 1.75, K = 0.0003. 

Simulating this system from the trivial initial condition leads to a limit cycle. In a continuation in K that 
limit cycle undergoes a subcritical period-doubling bifurcation. Now, we can do a two parameter continuation 
in (F, K) and draw the bifurcation diagram reported in Figure 10, 



5.4.1 The fold-flip point 

As can be seen in Figure 10, a fold-flip point is detected for (F, K) (1.7620, 0.2806 x 10~^). Since there is 
a Neimark-Sacker curve of the period doubled limit cycle rooted at the bifurcation point and the NS curve 
and the LPC curve lie on different sides of the PD curve, we are in the case represented in Figure 19 -(a), i.e. 
we have 020^11 < and ao2&ii < 0. Moreover, since the NS curve is supercritical, L^g should be negative. 
Numerically, we obtain that 611 = 562.2215, 020 = —0.5756, ao2 = — 0.1036e — 3, Lj^s = —178.8596. Hence, 
these results are in accordance with the theory. 
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Figure 10; Bifurcation diagram of a limit cycle in model (155). The blue curve is a limit point of cycles 
bifurcation, the green curve is a period-doubling curve (continue/dotted curve correspond to supercriti- 
cal/subcritical curves) and the magenta curve is a supercritical Neimark-Sacker bifurcation curve of the 
period doubled limit cycle. 



A Derivation of the normal forms 



A.l Notation 



Let M G M"^" be the monodromy matrix. In all codimension 2 cases all critical multipliers, i.e. all multipliers 
with modulus 1, have non-degenerate Jordan blocks. Let Mq be the critical Jordan structure, i.e. the block 
diagonal matrix consisting of the critical Jordan blocks, starting with the block of the trivial multiplier 1. 



Let /ife = e" 
defined as 



' (0 < 9k < Stt) be a critical multiplier with multiplicity m^.. The matrix G 







1 









1 

a-kj 



where ak is the Floquet exponent of the multiplier fik, with ak — iOk/T in the case of a positive real multiplier 
or a complex multiplier /i^ and Uk — ^ for /it = —1. The matrix Lq is the block diagonal matrix formed from 
the blocks Lk for which |/Zfc| = 1, starting with the block that corresponds with multiplier 1. The matrix Lq 
is the matrix Lq without the first row and the first column. 



A. 2 Bifurcations with 2 critical eigenvalues 
A.2.1 CPC 

At the CPC bifurcation the monodromy matrix has the critical Jordan structure 



Mo 



1 1 
1 
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A. 2 Bifurcations with 2 critical eigenvalues 



i.e. the multiplier /i = 1 is double non semi-simple. According to the proposed notation a — and thus 

^0 = ( n n ) ' -^0 0. 



We are in a situation where we can apply Theorem 2 from [23 . In particular, we have that 

where r plays the role of phase coordinate along the orbit. The polynomials p and P are T-periodic in r 
and at least quadratic in ^ such that 

^P{r, - ^Pir, OLoC = 0, ^Pir, + LIP{t, - ^Pir, - 

for all T and ^ G M. Putting Lq = we obtain 

ipir,O-0, iPir,O=0, 

i.e. the two polynomials are independent from r. So, by a Taylor expansion of the two polynomials the 
normal form becomes 

= 1 + e + p(o = 1 + e + aiC' + a'^e + • • • , 

= P{0 =be + r4^ + ... ■ 



Applying the chain rule gives 

dt 



+ c^^ + . . . )(1 + e + + "2^^ + ■ ■ ■ ) 



t3 



because of the cusp degeneracy condition. In order to obtain the proposed normal form (5), we do the 
substitution ^ i— )■ — ^ and find 

$ = 1 - ^ + "1^^ + "2^^ + . . . , 



dt 



with a2 = ^o^2- 



A.2.2 GPD 

At the GPD bifurcation the matrices described in Section A.l are 



1 

-1 








in = 0. 



We are in a case in which we can apply Theorem 3 from [23] . So we have the following 2T-periodic normal 
form (using the formula of Theorem 1 from [23] ) 



A Derivation of the normal forms 



67 



with polynomials p and P 2T-periodic in r and at least quadratic in ^ such that 

p{t + T, - p(t, -0, ^(t + T, ~0 = -^(^, 0- 

Putting io = in the first two formulas brings us in the same situation as of the previous case 

i.e. the two polynomials are independent of t. This makes it possible to rewrite the last two formulas as 

p(e) = p(-e), p{-o = -p{^i 

so polynomial p is even {p = 4>{i'^)) and polynomial P is odd {P = Therefore, taking the GPD 

degenerate condition into account, we can write down the first approximation of our normal form 

= l + 0(e') = l + aie'+a2e' + ---, 



Applying the chain rules gives 

because of the GPD degeneracy condition. So we obtain the normal form presented in (6), namely 



— = 1 + air + OiiC 
A. 3 Bifurcations with 3 critical eigenvalues 

A.3.1 CH 

In the CH case the Jordan block associated to the trivial multiplier is one-dimensional. We have 



'10 \ /O 

Mo = I e*"^ , Lo = | Lq 



ILO 

—iuj 



This puts us in a situation in which we can apply Theorem 1 from [53]. If we assume that ^ ^ Q, then it 
follows immediately from the results of Example III. 9 from j24j that the normal form is given by 
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A. 3 Bifurcations with 3 critical eigenvalues 



where the polynomials (/> and ij) are at least linear in their argument, cj) is real, while "0 takes values in C. 
We epand the polynomials up to the fifth order, namely 

dr 

- = l + ai|^|2 + a2|^|4 + .... 



dr 



such that the chain rule gives 



+ c'i\e + e'm^ + ...)(! + + «2iei' + . . . ) 



Here iujui + c' = ic since the Lyapunov coefficient of the Neimark-Sacker bifurcation is purely imaginary 
and iLoa2 + aic' + eJ — e. This gives us normal form (7). 

A.3.2 Rl 

At the Rl bifurcation the matrices described in Section A.l are 



'1 1 0\ /O 1 0> 

Mo = ( 1 1 , Lo = 1 I , Lo 
,0 1/ \0 0/ 



1 




We are in a case in which we can apply Theorem 2 from [23|. So we can define a T-periodic normal form 

J = l+6+p(r,0, ^ = h^ + P{r,0, 
where ^ = (^1,^2)- The polynomials p and P are T-periodic in r and at least quadratic in (^1,^2) such that 

If we write the polynomials in a Fourier expansion, namely 

00 00 



we obtain for any I € Z the following differential equations 

d 



d 



d^ 



C \ d(2) 



Putting our io into the equations and writing Pi{S,1t£,2) = (Pi {Cit£,2),Pi ('^11^2)) we can rewrite them as 
a set of differential equations in variable ^2 

"?2 J 4i 
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Since Pi{^i,^2), Pi '''\^i,^2) and P/^''(Ci,^2) are polynomials, if ^ ^ the only solution is the trivial one. So 
I equals zero and thus the polynomials are r independent. We obtain 



d 



0, 



The first two equations show that po and Pq^^ are independent from ^2, thus 

Po(6) = 0o(6), n^''(a) = eix(a). 

Integrating the last differential equation gives 

Po^'^(ei,6) = 6x(a) + ^(a). 

Now we can further simplify our normal form. In fact, since the homological operator associated with Lq 
has a two-dimensional null-space, we can make a change of variables such that polynomial Pg^^ vanishes (see 
[21] )■ So we can write 

where 0i and 02 are polynomials satisfying 0i(O) = 02(0) = = 0. 

Assembling all the information gives us the following normal form 

dr 

— = 1 + Ci + MCi) - 1 + 6 + a^i' + • • ■ , 

^ = i2Mii) + Mil) = ae, + 666 + • • • • 

Note that the polynomials 0o and 02 are at least quadratic in 61 while 0i is at least linear in its argument. 
In order to obtain the normal form presented in (8) we apply the chain rule which gives 



^ = 6(i + 6 + < + ...) 



and 



^ - («6' + 666 + ...)(i + 6 + "e? + ...) 

= ail + &C16 + . . . . 



A.3.3 R2 

At the R2 bifurcation the matrices described in Section A.l are 



'10 
Mo = ( -1 1 
.0 -1, 



'0 0\ 
Ln = I 1 



in = 



1 




^0 Oj 

We are in a case in which we can apply Theorem 3 from [23^ So we have a 2r-periodic normal form 

^ = l+p(r,6, §=^e + P(r,6, 
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A. 3 Bifurcations with 3 critical eigenvalues 



where ^ — (^1,^2)- The polynomials p and P are 2T-periodic in r and at least quadratic in their argument 
with 



(156) |:P(r,O-^P(r,C)iSe = 0, 

(157) p(t + T,0=p(t,-C), 



F(t + T,-0 = -P(t,0- 



Similar as in the Rl case (since the Lq matrix is the same) we obtain that all polynomials are independent 
from T, I has to be equal to zero and the first two polynomials p and P^^-* are independent from ^2- 
Since the polynomials obtained are independent from r, we can rewrite (157) as: 

p(0 = p(-0, P(-0 = -^(0, 

obtaining that polynomial p is even (p(^i) — ?!'o(Ci)) ^'^^ polynomials P(^i,^2) are odd (P'-^HCi) = 
Ci<^i(Ci) and P^^HCi; C2) = ^2^i(Ci) + Ci02(^i))- Now we can simplify our normal form by changing variables 
as discussed in the previous section. So we can write 

p(i)(a) = 0, p(^)(Ci,6) = 6'/'i(C?) + 602(€?)- 

Putting all information in the normal form equations gives the system 

dr 

— = 1 + 0o(Ci) = 1 + + . . . , 



^=6</'i(??)+ei02(C?)=ae? 



In order to obtain normal form (9) we apply the chain rule 



dt 



6 + aef 6 + . . . , 

«? + &C?6 + -..)(i + «^? + ---) 
+ + . . . . 



A.3.4 R3 

This is a simple case, since the Jordan block associated with the trivial multiplier is one-dimensional and -1 
is not a multiplier of the critical limit cycle. So we can write 



Mo 





Lo 



i^rj. 



We are in a case in which we can apply Theorem 1 from [23 . So we can define the following T-periodic 
normal form 

, , . dz - N 
— = 1+p{t,z), —^Loz + P{t,z), 

where z — (zi,zi). The polynomials p and P are T-periodic in r and at least quadratic in their argument 
such that 

|:P(t, z) - ^p(t, z)L*z = 0, ^P(t, z) + L*PiT, z) - ^^P{t, z)L*z = 0. 



A Derivation of the normal forms 



71 



We apply the results obtained in Example III. 9 from [53] with wT/27r = 1/3 to obtain 
at 

dzi l27r ^ /I |2 3 -i2TTT/T\ I -2 i2-!TT/T^ i\ |2 -3 i2-KT IT\ 

■ dr 61 

Defining a new variable ^ — e~^^^ zi, this system can be rewritten as 



with polynomials (/)o and (pi at least linear in their arguments, while 02(0) 7^ 0. Notice that this system is 
autonomous and equivariant under the rotations of angle 27r/3. Expanding the polynomials gives 



dr 

— = l + ai|e|'+a2e^ + a3f + 



dr 



and thus 



dt 



(fee' + cCleP + ...)(!+ + «2e' + «3f + . . . ) 



- fef + ce|^P + ..., 
so normal form (10) is obtained. 

A.3.5 R4 

As in the previous case the Jordan block associated with the trivial multiplier is one-dimensional. The 
matrices in Section A.l are 



^10 
Mo = ( e't 

,0 e-'t 



^0 



.0 



* 2T ^ 
~*2T 



We can apply Theorem 1 from [53] and define a T-periodic normal form 

dr dz ~ 

— = l+p{T,z), — = Loz + P{t,z}, 

at dr 

where z — (zi,zi). The polynomials p and P are T-periodic in r and at least quadratic in their argument 
such that 

^p{t, z) - ^p{T, z)L*z = 0, -f P(r, z) + L*P{t, z) - ^P{t, z)L*z = 0. 
dr dz dr dz 

Again, we use Example III. 9 from [24] with ujT/2tt = 1/4 and obtain 



^ = l + Po(kiP,^1e''^^/^,z4e-2-/^), 

^ = ^^l + ^lQo(|^l|^4e-^''^^/^) + ^1e''-^/^gl(|.l|^z1e^''^^/^). 
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A. 3 Bifurcations with 3 critical eigenvalues 



Defining a new variable ^ = e zi, the system can be rewritten as 

i + 0o(iep,f,a, 



'dt 

dr 



with polynomials 00 and 01 at least linear in their arguments, while 02(0) ^ 0. Notice that this system is 
autonomous and equivariant under the rotations of angle 7r/2. Expanding the polynomials gives 



dr 
dt 

dr 



i + «iiei' + «2e' + a3e' + 

ceieP + rff + . . . , 



We still have to transform the second equation, namely 



dt 



(c^l^P + + . . . )(1 + ail^p + a2^^ + ag^"^ + . . . ) 

cciei' + df + ..., 



so normal form (11) is obtained. 
A. 3. 6 LPPD 

At the LPPD bifurcation the matrices from Section A.l are 



Mo 










We are in a case in which we can apply Theorem 3 from 23 . So we can define a 2T-periodic normal 
form 

^-l + ei+p(r,0, ^ = Lo^ + P(t,0, 

where ^ — (^17^2)- The polynomials p and P are 2T-periodic in r and at least quadratic in their argument 
such that 



(158) 

(159) 
(160) 

By putting Lq into (158), we obtain 

d 



P(i) (r + T, 6,-6) = P^'^ (r, a , 6), P^'^ ir + T, ^i, -6) = -P^'^ (r, 6 , 6). 



ar ar dr 



0, 



so our polynomials are independent from r. Then, using (159) and (160), there holds 

P(ei,6) =p(a,-6), P^'\^i,-^2) = p(^HCi,6), P^'\^i,-^2) = -p('Hei,6), 

so the polynomials are of the following form 

P = Xl(ei)+X2(el)(l + X3(ei)), 

P^^^ = 6^1(^1) + ^2^2{il){l + </^3(ei)), 
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with xi s-iid V'l s-t least quadratic in their argument and all the other polynomials at least linear in their 
argument. 

Assembling all the information gives the following system 



dr 
'dt 

d^ 
dT 

dr 



= l+ei+Xi(ei)+X2(e|)(l + X3(ei)) 

= 1 + a + "20^? + "02^1 + asofi + a'i2^i£.i 

= 6<^i(6) + 6^2(el)(i + '^3(a)) 

= &iia6 + fo$>ie?6 + bosCi + ■■■■ 



Applying the chain rule gives 



da 

dt 



with ago = a!jo + 4o'«i2 



= 1^20^1 + "02C2 + ^SO?! + "20^1 + '^12^1^2 + '^02^1^2 + 

= a2oCi + ao2^2 + asoCi + ai2a^2 + ■ • • : 
and 



''12 



-'02! 



with 621 = 621 + ^'11 ■ N^ow, we do the substitution a ^a which gives 

dr 



= 1 - ^1 + "20^1 + ao2'?| - "30^1 - "126^2 

da 

t 



'^a _ / c2 ' c2 I (-3 I c c2 

— —020?! — ao2S2 + 030?! + ai2?l?2 



By putting a^o = -0:30,012 = -Q;i2,a2o = -030,002 = -002,^11 = -^'n, we obtain (12). 



B Poincare maps of the periodic normal forms 

B.l Bifurcations with 2 critical eigenvalues 
B.l.l CPC 

If we reparametrize time, we obtain the system 

1, 

ce 

2 h • • • . 

1-^+01^2+02^^ 

Note that we can define a Poincare map out of system ( 5 ) evaluating the solution of this system starting 
from a point (0,77) at time t = T. Since the two equations are uncoupled we can consider only the second 
one. Making one Picard iteration [29], we can construct an approximation of this map as follows 

Co = a = Co + / 

Jo 1 - Co + aiCo + "2C0 
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B.l Bifurcations with 2 critical eigenvalues 



and expanding in a Taylor series we obtain 

(161) r]>-^f] + cTrf + 0{ri^). 

Further iterations do not change this expansion and so (161) is the Poincare map of system (5). Note that 
this Poincare map is similar to the one for the cusp point of fixed points [5^. Therefore, we can conclude 
that the behavior of the system in the neighborhood of the bifurcation is the same. In particular, referring 
always to [35], we can draw the bifurcation diagram of the Poincare map and obtain Figure 11 On the two 
drawn curves, with label Ti and T2, two limit cycles collide and disappear. The output given by MatCont 
is the normal form coefficient c. 




(a) c < (b) 



Figure 11: Bifurcation diagram of the cusp bifurcation of the fixed point normal form. 



B.1.2 GPD 

First, we have to calculate the second iterate of the normal form of the Generalized flip bifurcation of fixed 
points: 

(162) 5(„) = _^ + _Lg^;5 ^ g2(^) _Lg^5^ 

Reparametrizing the time of (6) and doing Taylor expansion up to the fifth order gives 




where the dots are O(C^) terms 2T-periodic in r. Doing one Picard iteration of the second equation up to 
2T we obtain the Poincare map of our normal form 

e ^ e + "^Tee + ■■■ 

and see that it is the same map as the second iterate of the generalized period-doubling normal form of fixed 
points. In particular, since the coefficient of the fifth order term of the second iterate of the normal form of 
fixed points has opposite sign than the one of limit cycles, we can conclude that the behavior of the system 
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at the bifurcation is the same but with opposite sign of the normal form coefheient. Moreover, since T is 
nonzero, the non-degeneracy condition is e ^ 0. In fact, if e > we obtain the bifurcation diagram reported 
in Figure 12 -(a), in which the hmit point bifurcation of the period doubled limit cycles (T*^^)) is tangent to 
the supercritical period-doubling branch (the one in which the normal form coefficient is positive, so labeled 
as F^^), while if e is negative we are in the opposite situation, reported in Figure 12 -(b). The output given 
by MatCont is the normal form coefficient e. 




(a) e > (b) e < 



Figure 12: Bifurcation diagram of the degenerate period-doubling point bifurcation of the fixed point normal 
form. 



Remark that the normal form coefficient of the period-doubling bifurcation of limit cycles, computed 
through periodic normalization as done in [30], has opposite sign of the one for fixed points. Also in this 
case, the normal form coefficient of the GPD bifurcation of limit cycles has opposite sign of the one for fixed 
points. 



B.2 Bifurcations with 3 critical eigenvalues 

B.2.1 CH 

In this section we will show how the periodic normal form (7) is related to the normal form of the Chenciner 
bifurcation of fixed points, and how the non-degeneracy condition of the Chenciner bifurcation of limit cycles 
is related to the Chenciner bifurcation of maps. First note that, if we scale the time, we can rewrite system 
(7) as 

— = zw^ + i(c - q;iw)^|^P + (e - «(aic - a\LU + Q!2a;))^|^|'' + . . . . 
■ dr 

We need to make a change of variables in order to obtain a quasi-identity flow. Introducing the new 
complex variable z — e~'"'^^, the second equation can be rewritten as 

dz 

— — i{c — aiuj)z'^z + (e — i(aic — aiuj + a2w)) z^z"^ + . . . . 
dr 

Doing two Picard iterations up to time T, we obtain the rotating Poincare map of the system 

/ r^T 1 \ - 

C ^ C + iT{c - aia;)^2^ -|- T f e - — -|- aicTw - -aiTw^ -|- i (a?a; - aiC - aaw) j i^i^ + ... . 
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B.2 Bifurcations with 3 critical eigenvalues 



Note that this system has the same Poincare map as the Chenciner bifurcation in the fixed point case [IH] , so 
we expect the same bifurcation scenario on the Poincare map of the system. Notice that the real part of the 
first Lyapunov coefficient is 0, since the Neimark-Sacker bifurcation is degenerate. The sign of the second 
Lyapunov coefficient L2 (as defined on page 420 of [29]) determines the bifurcation scenario. However, from 
(7) we can derive that 3?(e) < corresponds with a stable critical limit cycle and 3i{e) > with an unstable 
critical limit cycle. Therefore, the case 3?(e) < corresponds with the case L2 < and ^R{e) > corresponds 
with L2 > 0. So 3fi(e) and the second Lyapunov coefficient L2 as defined in [29] have the same sign and 
vanish at the same time. Since both coefficients have the same effect and L2 requires more computations, we 
compute 3?(e) to determine the bifurcation scenario, and in this paper we will call 3?(e) the second Lyapunov 
coefficient. The bifurcation diagram in the neighborhood of this codim 2 point can be found by looking at 
the sign of this value. When 5R(e) < the outer invariant curve is stable and the limit point of tori curve 
is tangent to the subcritical Neimark-Sacker branch, as shown in Figure 13 -(a). When 3?(e) > the outer 
invariant curve is unstable and the limit point of tori curve is tangent to the supercritical Neimark-Sacker 
branch, see Figure 13-(b). The output given by MatCont is 3fJ(e). 




(a) K{e) < (b) 5R(e) > 



Figure 13: Bifurcation diagram of the generalized Neimark-Sacker bifurcation of the fixed point normal form. 



B.2. 2 Rl 

In this section we will show how the periodic normal form (8) is related to the normal form for the 1:1 
resonance of fixed points, which allows us to formulate the non-degeneracy condition. First note that, if we 
scale the time, we can rewrite system (8) as 

1, 

i-Ci + < ' 



'dr _ 

da _ 
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Doing Picard iterations up to time T, it's possible to show that the generated map is the same, up to second 
order terms, of the one generated by Picard iterations up to time 1 of the truncated system 

= r(6 + 2^6), 

This last system is topological equivalent to the one-time shift ODE which represents the first iteration of 
the normal form of 1:1 resonance of fixed points, i.e. the system 

Cl = C2, 

C2 =aiC?+6iCiC2, 

since we can transform one in the other using the following change of variables 

The bifurcation phenomena in the 1:1 resonance of cycles are the same of those which appear in the 1:1 
resonance of fixed points. Notice that 

ai — T^a, bi = Tb, 

so the non-degeneracy conditions are 

a 7^ 0, b^O, 

and the cases depend on the sign of the product of a and b. In particular, as shown in figure 14, if the two 
coefficients have different sign the Neiniark-S acker bifurcation (labeled H) is supercritical (with negative 
normal form coefficient), while in the other case it is subcritical. The output given by MatCont is the 
product of the coefhcients a and b. 




(a) a> 0, ab <0 (b) a > 0, afe > 



Figure 14: Bifurcation diagram of the 1:1 resonance bifurcation of the fixed point normal form. The other 
two cases in which a < can be obtained by a reflection around the origin of the state portraits and a 
horizontal flip of the bifurcation diagrams. 
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B.2 Bifurcations with 3 critical eigenvalues 



B.2.3 R2 

In this section we will show how the periodic normal form (9) is related with the normal form for the 1:2 
resonance of fixed points, which allows us to formulate the non-degeneracy condition. First note that, if we 
scale the time, we can rewrite system (9 ) as 

1, 

l + < 

1 + <2 + • • • • 

Doing Picard iterations up to time 2T, it's possible to show that the generated map is the same, up to third 
order terms, of the one generated by Picard iterations up to time 1 of the truncated system 

This last system is topological equivalent to the one-time shift ODE which represents the second iteration 
of the normal form of 1:2 resonance of fixed points, i.e. the system 

Ci = C2, 

C2 = aiC?+6iCfC2, 

since we can transform one in the other using the following change of variables 

C2 = 

The bifurcation phenomena in the 1:2 resonance of cycles are the same that appear in the corresponding 
bifurcation of fixed points. Notice that in this case 

ai = AT'^a, 61 = 2Tb. 

so the non-degeneracy conditions are 

a ^ 0, by^O. 

We have four different unfoldings as possible bifurcation diagrams, determined by the signs of the coefficients. 
The ones with negative b are reported in Figure 15 The other two cases can be obtained by reversing the 
arrows of the phase portraits and making a vertical flip both of the state portraits and of the bifurcation 
diagrams. The primary Ncimark-S acker bifurcation (labeled iJ*^^') is supercritical (with negative normal form 
coefficient) if our coefficient b is negative, subcritical otherwise. Moreover if a < a secondary Neimark- 
Sacker bifurcation (H^^^) is rooted at the 1:2 resonance point with opposite criticality of the primary one. 
The output given by MatCont is (o, b). 

B.2.4 R3 

In this section we will show how the periodic normal form (10) is related with the normal form for the 1:3 
resonance of fixed points, which allows us to formulate the non-degeneracy condition. First note that, if we 
scale the time, we can rewrite system (10) as 



dT 

, dr 
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(a) 6 < 0, a > (b) 6 < 0, a < 



Figure 15: Bifurcation diagram of the 1:2 resonance bifurcation of the fixed point normal form. The other 
two possible cases in which b > can be obtained by reversing time and making a vertical flip both of the 
state portraits and of the bifurcation diagrams. 



Doing Picard iterations up to time 3T, it's possible to show that the generated map is the same, up to cubic 
terms, of the one generated by Picard iterations up to time 1 of the truncated system 

i = 3Tb£^ + 3Tc4\^\^. 

This system is topological equivalent to the one-time shift ODE which represents the third iteration of the 
normal form of 1:3 resonance of fixed points, i.e. the system 

The bifurcation phenomena in the 1:3 resonance of cycles are the same of the ones which appear in the 1:3 
resonance of fixed points, if the non-degeneracy conditions are satisfied, i.e. 

b ^ 0, 3?(c) ^ 0. 

As can be seen in Figure 16, if 3fi(c) < the Neimark-S acker bifurcation (labeled N) is supercritical (with 
negative normal form coefficient), while in the other case it is subcritical. The output given by MatCont is 
(6,5i(c)). 



B.2.5 R4 



In this section we will show how the periodic normal form (11) is related with the normal form for the 1:4 
resonance of fixed points, which allows us to formulate the non-degeneracy condition. First note that, if we 
scale the time, we can rewrite system (11) as 

|=cC|Cp-fdf + .... 

Doing Picard iterations up to time 4T, it's possible to show that the generated map is the same, up to cubic 
terms, of the one generated by Picard iterations up to time 1 of the truncated system 
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B.2 Bifurcations with 3 critical eigenvalues 




(a) 5R(c) < (b) 5R(c) > 



Figure 16: Bifurcation diagram of the 1:3 resonance bifurcation of the fixed point normal form. 

This system is topologically equivalent to the one-time shift ODE which represents the fourth iteration of 
the normal form of 1:4 resonance of fixed points, i.e. the system 

C = ciCICI' + rfiC'- 

The bifurcation phenomena in the 1:4 resonance of cycles are the ones which appear in the corresponding 
bifurcation of fixed points, if the non-degeneracy condition is satisfied, i.e. 

After defining 




we can determine which bifurcations occur by looking at the place of A in the Gauss plane (see Figure 17). 
We here report the different possible bifurcation diagrams. A complete investigation can be found in the 
literature P?!!^- 

Many topologically different bifurcation diagrams can be found near a 1:4 resonance point. The analysis, 
if one excludes higher codimension situations, can be reduced to 22 different cases, which, as mentioned 
before, depend on the value of A. First of all, analyzing the normal form, one can divide the Gauss plane 
into two big regions: in the semiplane < the primary Neimark-Sacker bifurcation is supercritical, in 

the semiplane ^{A) > it is subcritical. What happens in the semiplane > can therefore be obtained 

by inverting the direction of the vector fields. We can further reduce the analysis to the third quadrant of 
the Gauss plane, since the 12 possible cases are topologically equivalent paired through the transformation 
C I— > C. The different regions are reported in Figure 17, in which only some curves (the continuous lines) are 
known analytically, the dashed curves are computed numerically. 

Figure 18 shows the possible bifurcation diagrams with the sketches of the phase portraits for the Poincare 
maps in the case that ^{A) < 0. Many local and global bifurcations are involved in the different scenarios. 
We use the following notation, consistent with the rest of the text: 



N: Neimark-Sacker bifurcation. In regions VII and VIII we also have a Neimark-Sacker bifurcation of the 
period-4 limit cycle. 
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-4-3-2-1 1 2 3 4 



Figure 17: Partitioning of the A plane into topologically different regions. 



T: Fold bifurcation of the period-4 limit cycles. There are three possibilities. Superscript in, on or out 
means that the bifurcation happens inside, on or outside the invariant curve. 

H: Homoclinic connection of the period-4 saddle limit cycle. Superscript S means that the born invariant 
curve is smaller than the limit cycle (a square looking homoclinic connection), C that it is bigger (a 
clover looking homoclinic connection), and L means that the born invariant curve is around the period- 
4 limit cycle; subscript + (-) means that the saddle quantity is positive (negative), so the horning 
invariant curve is repelling (attracting). 

F: Fold bfirucation of the tori. 



The output given by MatCont is {A,d). 




Figure 18: Bifurcation diagrams locally to the 1:4 resonance bifurcation in the different regions of figure 17 
The cases in which Re{A) > can be obtained with the transformation t — > —t, (3 — ?> —(3. 
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B.2.6 LPPD 

In this section we will show how the periodic normal form (12) is related to the normal form for the fold-flip 
of fixed points, which allows us to formulate the non-degeneracy condition. As for the R2 case we show the 
topological equivalence of the 2r-shift map of our system with the 1-shift map of the approximating vector 
field of the fold-flip normal form for fixed points. First of all we can prove, by scaling the time and using 
Picard iterations, that the Poincare map of our normal form is the same, up to cubic terms, as the one-shift 
map of the truncated system 

= 2Ta2o(l + 2^002^2 + 2^(030 + 020)^? + 2T(ai2 + ao2)^i^l 
\C2 = 2T6na6 + 2T(62i + 5u)e?6 + ^Tbos^l 

This system is topologically equivalent with the one-time shift ODE which represents the second iteration 
of the normal form of the fold-flip bifurcation of fixed points, i.e. the system 



(163) 



Ci = aid + bid + (ci - al)Cf + (di - ai6i -f 6i)CiC2', 
C2 = -CiC2 + ^(ai-l)CiC2 + J^^iCa', 



since this second system can be obtained (neglecting higher order terms) from the first one using the trans- 
formation 

fCi - -2hiT^i - 2T (&n + 621 + a2obiiT + bj^T) ^ - 2T(&03 + ao^biiT)^, 
\C2 = 26nrC2. 

This transformation should be invertible, so one non-degeneracy condition is involved, namely 

611 0. 

If this condition is satisfied, then the system can be put in the form (163), where the constants are defined 
as 

fl20 ao2 _ ^20 + ^30 + 2a2QT 



ai = --; — , bi = — , ci 

dl : 



611' 611' 2bl,T 

2a2ofoo3 + 3ao26ii + ai26ii + 25o3feii + 2ao262i + 2002020^11^ + 60026?!^ 



2bl,T 

and from those values we can understand which types of bifurcation the system has. In particular (see 
[291 133| for more details) three more non-degeneracy conditions are involved 

• if ^20 7^ there are two limit cycles that collide and disappear 

• if ^02 7^ a period doubled limit cycle born in this point. 

Moreover if ^ < a torus bifurcation occurs on the period doubled orbit, with Lyapunov coefficient 

Lns = -2a2o&03 - 3002030^11 + a2o(ai2&ii + 26o3&ii + 2ao2&2i) 

and so, in order to avoid degeneracy, we also assume Lns 7^ 0. 

In Figure 19 the four possible scenarios are reported depending on the sign of the normal form coefficients. 
The output given by MatCont is (611, 020, ao2, ^Afs)- 



84 





(D 


@ 






( ® 




k 






-y® 







(2) 


@ 






I ® 




(o P 






. ® 
p_ 





L*<s> 



Us< 



© 





(a) 020611 < 0, 002611 < 



(b) 020611 < 0, 002611 > 

(D 



® 





® 







@ \ 




NS, 


® 


® ^) 




^'-^^ 








Lns 


> 











® 


Q) 












p_ 


Lns 


< 




^1) 






® 


(D J 














(c) 020611 > 0, 002611 < (d) 020611 > 0, 002611 > 

Figure 19: Bifurcation diagrams of a fold-flip bifurcation of the fixed point normal form. 

C Kernels of some differential-difference operators 

In Section 4 wc used the orthogonality with respect to the following inner product: if Ci)C2 € C°([0, 1],C") 
and ?7i, ?72 € C", then 

If this inner product vanishes, then we say that the corresponding vectors are orthogonal and write 







■ C2 " 






. ^2 _ 



In Section 4 we used some propositions, from which we will give the proof in this appendix. 
Proposition C.l. Consider two differential- difference operators 

.^1,2 :Ci([0,l],R") ^C°([0,1],R") xR", 
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with 



C-TAit)C 

C(o)-C(i) 



,M0 



M0 = 

//C e Ci([0, then C e Ker(^i) if and only if 
and C € Ker(^2) */ and only if 



C(o)-C(i) 



c 

C(o) 



±0i(C^([O,l],M")). 



Proof. We will focus on the first assertion. If ^ is in the kernel of (/>i , then ( — TA{t)( = and C(0) — C(l) = 0- 
For all g e C^([0, 1],M") we have 

g'^mm Tg'^{t)A{t)at)dt = O 

^ 9''mt)\i - lo g'^mm - Tg'^{t)Amt)dt = o 

^ 5^(l)C(l) - .9^(0)C(0) - /o (gW + TA^{t)g{t)rat)dt = 

(5(0) - gii)rm - j^im + TA'^{t)g{t)rcm = o 



g + TA^{t)g 
9(0) -5(1) 



C(o) 



0. 



Conversely, assume that ( 



" c ■ 




'g + TA'^{t)g] 


C(0). 




. 5(0) -5(1) J ^ 



) = for all g e C\[0, 1],K"). Then, 



/o + TA^{t)g{t))dt + C^(0)(.g(0) - .9(1)) = 

=> C{l)g{l) (0)5(0) - /„\C(t) - TA{t)mfg{t)dt + C^(0)(5(0) - 5(1)) = 

^ -(c(o) - c(i))^5(i) - Sim - TAmWgm = o. 

If C{t) - TA{t)Cit) ^ 0, then there exists a g{t) with 5(1) = such that 



m-TA{t)at)fg{t)dt^o. 



This is impossible, so ({t) - TA{t)Cit) = 0. Hence (C(0) - C(l))'^.g(l) = for all g; and thus there must hold 
that C(0) — C(l) = 0. From both observations it follows that ( G Ker(^i). 

The proof of the second assertion is similar. □ 

Proposition C.2. Consider ^1,2 : C^{[0, 1],K") ^ C°([0, 1],M") x M", where 



C - TA{t)C 
C(0)+C(1) 



,02(0 



<^i(C) = 

//C e Ci([0, 1],R"), then C S Ker(0i) if and only if 

^1^3^ ]±02(C^([O,1],M")) 

and C, S Ker(02) if and only if 



c+TA^m 

C(o) + C(i) 



c 

C(o) 



Proof. The proof is similar to the proof of Proposition C.l, 



□ 
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Proposition C.3. Consider (t>i,2 ■ C^{[Q, 1], C") C°([0, 1], C") x C", where 



i(C) 



C - TA{t)c + lec 
C(o)-C(i) 



,02(0 



^ + TA^{t)c + iec 
C(o)-C(i) 



//C e Ci([0, 1],C"), i/ien C e Ker(0i) i/ and onZy i/ 

^^^Q^ ±</.2(Ci([0,l],C")), 

anrf C € Ker(02) */ on/y i/ 



c 

C(o) 



±,^i(Ci([0,l],C")). 



Proo/. If C is in the kernel of (t>i, then C - TA{t)C + i6'C = and C(0) - C(l) = 0. For all g £ C\[0, 1], C") 
we have 

g^{t)cm Tg^marndt + leg^mt) = o 
^ g'^mm ~ lo g'^mm - j^' Tg^mit^m + zog^mt) = o 

^ 5"(1)C(1) - 5"(0)C(0) - J^m + TA^{t)g{t) + ieg{t)rC{t)dt = 
-{g{0) - 5(1))"C(0) - /o + TA-^mt) + ieg{t))^at)dt = 



g + TA'^{t)g + i0g 
.9(0) -.9(1) 



c 

C(o) 



0. 



The proofs of the reverse implication and the second assertion are similar. 
Proposition C.4. Consider ^1,2 : C^HO, 1],C") C°([0, 1],C") x C", where 



C-TA{tK 
C(0) - e-«C(l) 



,M0 = 



MO = 

//C e Ci([0, 1], C"), then C G Ker(0i) if and only if 

^^^Q^ ]±02(C^([O,1],C")) 

anrf C € Ker(^2) if o-nd only if 



C(0) - e-«C(l) 



c 

C(o) 



±<^i(Ci([0,l],C")). 



□ 



Proo/. If C is in the kernel of (^i, then C - TA{t)C = and C(0) - e-'^C(l) = 0. For all g e C^{[0, 1],C") we 
have 

/o g'^WCm - Tg^{t)Amt)dt = 
^ 5"(t)C(t)lJ - /o ff"(OC(Odi - /o Tg^mmt)dt = 
^ 5"(1)C(1) - ff"(0)C(0) - /o (5(t) + TA'^{t)g{t))^at)dt = 
-(5(0) - e-^Og{l)mO) - /o (.9(i) + TAT(f)5(t))H^(f)df = 
g + TA^{t)g C 
5(0)-e-^(l) ' C(0) 



0. 



The proofs of the reverse implication and the second assertion are similar. 



□ 



C Kernels of some differential-difference operators 



87 



Proposition C.5. Consider two differential- difference operators 0i,2 : C^([0, 1],]R") C°([0, 1],R") x R", 
where 



MO 



C - TA{t)C 

C(o)-C(i) 



c 




. C(o) . 





01 (0 



h + TA^{t)h 
h{0) - h{l) 



2(C) 



c+TA^m 

C(o)-C(i) 



9 




■ h ' 











IfC&C\[0,l],W), then 



if and only if 



for all h e C^([0, Furthermore 



if and only if 



for allheC^{[0,l\,W). 

Proof Wc focus on the first assertion. Suppose that C{t) - TA{t)C{t) = g{t) and C(0) - C(l) = 0. For all 
h e Ci([0,l],R") we have 

/o h^{t)t{t)dt - Th^{t)A{t)C{t)dt = h^{t)g{t)dt 
h^{t)m\h - lo h^{mt)dt Th^{t)A{t)at)dt = h^{t)g{t)dt 
^ /iT(l)C(l) - /iT(0)C(0) - jl{h{i) + TA^{t)h{t)Yat)dt = /o h^{t)g{t)dt 
Jo C{t){h{t) + TA^{t)h{t))dt + C'^(0)(/i(0) - h{l)) = - Jl g'^{t)h{t)dt 
h + TA^{t)h" 



c 




■ h - TA{t)h ' 




9 




■ h ' 


. C(o) . 




h{0) - /i(l) 













c 

C(o) 



/i(0) - /i(l) 

The proofs of the reverse implication and the second assertion are similar. 



□ 



Proposition C.6. Consider two differential- difference operators 4>i^2 ■ Ci([0, 1],E") C°([0,1],R") x M", 

where 



MO 

//C eCi([0,l],R"), then 
if and only if 
V/ieCi([0, 1],IR"). Furthermore 
if and only if 
V/ieCi([0,l],K"). 



C - TA{t)C 

C(o) + C(i) 



, 02(0 



C + TA-^itX 

C(o) + C(i) 



c 




. C(o) . 


7 



01 (0 = 



h + TA^{t)h 
h{0) + h{l) 



9 




■ h ' 





7 






02(C) 



c 




■ h - TA{t)h ' 




9 




■ h ' 


. C(o) . 




_ /i(0) + /i(l) . 
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Proof. Suppose that C(0 - TA{t)Cit) = g{t) and C(0) + C(l) = 0. For all h G C^{[Q, 1],]R") we have 



K^{t)at)dt - Th^{t)A{t)at)dt = h^{t)g{t)dt 

/iT(l)C(l) - /i^(0)C(0) - !l{h{t) + TA^it)hit)rcit)dt = h^{t)g{t)dt 
J^Cit)ihit) + TA^it)hit))dt + Cmh{0) + hil)) = -J^g'^{t)h{t)dt 



c 

C(o) 



h + TA^{t)h 
h{0) + h{l) 



The proofs of the reverse implication and the second assertion are similar. 



□ 
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